# Stochastic Unit Commitment in Low-Inertia Grids

Matthieu Paturet, Uros Markovic, Stefanos Delikaraoglou,
Evangelos Vrettos, Petros Aristidou, and Gabriela Hug

###### Abstract

In this paper, the Unit Commitment (UC) problem in a power network with low
levels of rotational inertia is studied. Frequency-related constraints, namely
the limitation on Rate-of-Change-of-Frequency (RoCoF), frequency nadir and
steady-state frequency error, are derived from a uniform system frequency
response model and included into a stochastic UC that accounts for wind power
and equipment contingency uncertainties using a scenario-tree approach. In
contrast to the linear RoCoF and steady-state frequency error constraints, the
nadir constraint is highly nonlinear. To preserve the mixed-integer linear
formulation of the stochastic UC model, we propose a computationally efficient
approach that allows to recast the nadir constraint by introducing appropriate
bounds on relevant decision variables of the UC model. For medium-sized
networks, this method is shown to be computationally more efficient than a
piece-wise linearization method adapted from the literature. Simulation
results for a modified IEEE RTS-96 system revealed that the inclusion of
inertia-related constraints significantly influences the UC decisions and
increases total costs, as more synchronous machines are forced to be online to
provide inertial response.

###### Index Terms:

Unit commitment, low-inertia grid, frequency constraints, wind uncertainty,
voltage source converter.

## Nomenclature

The main notation used in this paper is introduced below. Additional symbols
are defined in the paper where needed. All symbols are augmented by index $t$
when referring to different time periods.

###  -A Sets and Indices

  * $\ell\in\mathcal{L}$

Set of transmission lines.

  * $\xi\in\mathcal{E}$

Set of scenarios $\xi=\\{c,\omega\\}$ including generation outages ($c$) and
wind power uncertainty ($\omega$).

  * $i\in\mathcal{I}$

Set of conventional generation units.

  * $j\in\mathcal{J}$

Set of converter-based (i.e., wind) generation units.

  * $n\in\mathcal{N}$

Set of nodes.

  * $\mathcal{I}_{n}$

Set of conventional generation units located at bus $n$.

  * $\mathcal{J}_{n}$

Set of converter-based units located at bus $n$.

###  -B Decision variables

  * $\hat{\delta}_{n}$

Day-ahead voltage angle at node $n$ $[\mathrm{rad}]$.

  * $\tilde{\delta}_{n\xi}$

Real-time voltage angle at node $n$ in scenario $\xi$ $[\mathrm{rad}]$.

  * $F_{\xi t}$

Global fraction of total power generated by high-pressure turbines in scenario
$\xi$ $[\mathrm{p.u.}]$.

  * $k_{i\xi t}$

Scaled power gain factor of conventional unit $i$ $\xi$ $[\mathrm{MW}]$.

  * $l_{n\xi t}^{\mathrm{shed}}$

Shedding of load at node $n$ in scenario $\xi$ $[\mathrm{MW}]$.

  * $M_{\xi t}$

Global system inertia in scenario $\xi$ $[\mathrm{p.u.}]$.

  * $p_{it}$

Day-ahead dispatch of conventional unit $i$ $[\mathrm{MW}]$.

  * $r_{i\xi t}^{+/-}$

Up-/Downward reserve deployment of unit $i$ in scenario $\xi$ $[\mathrm{MW}]$.

  * $R_{\xi t}$

Global system droop factor in scenario $\xi$ $[\mathrm{p.u.}]$.

  * $u_{it}$

Commitment variable of conventional unit $i$.

  * $w_{jt}$

Day-ahead dispatch of wind power unit $j$ $[\mathrm{MW}]$.

  * $w_{j\xi t}^{\mathrm{spill}}$

Wind spillage of unit $j$ in scenario $\xi$ $[\mathrm{MW}]$.

  * $y_{it}$

Start-up variable of conventional unit $i$ .

  * $z_{it}$

Shut-down variable of conventional unit $i$ .

###  -C Parameters

  * $\alpha_{i\xi t}$

Outage parameter of conventional unit $i$ in scenario $\xi$.

  * $\Delta P_{\xi t}$

Size of power outage in scenario $\xi$ $[\mathrm{p.u.}]$.

  * $\pi_{\xi}$

Probability of occurrence of scenario $\xi$.

  * $B_{nm}$

Susceptance of transmission line $(n,m)$ $[\mathrm{S}]$.

  * $C_{i}$

Day-ahead price offer of unit $i$ $[\mathrm{\$/MWh}]$.

  * $C_{i}^{\mathrm{SU/SD}}$

Start-up/Shut-down price offer of unit $i$ $[\mathrm{\$}]$.

  * $C_{i}^{\mathrm{+/-}}$

Up-/Downward reserve price offer of unit $i$ $[\mathrm{\$/MWh}]$.

  * $C^{\mathrm{sh}}$

Value of lost load $[\mathrm{\$/MWh}]$.

  * $D_{nt}$

Demand at node $n$ $[\mathrm{MW}]$.

  * $\overline{f}_{nm}$

Capacity of transmission line $(n,m)$ $[\mathrm{MW}]$.

  * $\overline{P_{i}}/\underline{P_{i}}$

Active power limits of conventional unit $i$ $[\mathrm{MW}]$.

  * $R_{i}^{\mathrm{+/-}}$

Up-/Downward reserve capacity of unit $i$ $[\mathrm{MW}]$.

  * $R^{\mathrm{U}/\mathrm{D}}_{i}$

Ramp up/down limits of conventional unit $i$ $[\mathrm{MW/h}]$.

  * $W_{j\xi t}^{*}$

Wind power realization of unit $j$ in scenario $\xi$ $[\mathrm{MW}]$.

  * $W_{d}$

Total capacity of droop control units $[\mathrm{MW}]$.

  * $W_{v}$

Total capacity of VSM units $[\mathrm{MW}]$.

##  I Introduction

With increasing penetration of renewable energy sources, system operators face
new challenges in order to ensure power grid stability. One of these
challenges is frequency stability due to a loss of generation or a large
variation of load. In traditional power systems, synchronous generators (e.g.,
hydro or steam turbines) provide rotational inertia through stored kinetic
energy in their rotating mass (turbine system and rotor). This energy is
important to stabilize the system as it ensures slower frequency dynamics and
reduces the Rate of Change of Frequency (RoCoF) in case of a generation-demand
imbalance [1]. In the future, with more generation coming from wind and solar
power, the ability of the system to maintain the frequency within the
acceptable range is diminished. Indeed, photovoltaic systems are connected to
the grid through inverters, which do not exhibit rotational inertia. Even in
the case of inverter-interfaced wind generators, the inverter electrically
decouples the rotor’s rotational inertia from the system [2].

Transmission System Operators (TSOs) around the world are concerned with the
stability issues associated with large penetration of renewable energy in
their systems. In the United States, the Electric Reliability Council of Texas
(ERCOT) has studied the effect of low inertia on the security and reliability
of the grid [3]. The Irish TSO, EirGrid, is designing ancillary services to
remunerate providers of rotational or synthetic inertia [4]. Furthermore,
EirGrid currently imposes limits on the maximum instantaneous penetration of
variable RES with respect to the total load demand at any point in time.

In systems with low rotational inertia, TSOs must impose minimum inertia
requirements in order to secure frequency stability and avoid system collapse
in case of a severe fault or a sudden mismatch between generation and demand.
With such new requirements, the traditional Unit Commitment (UC) problem,
i.e., the day-ahead scheduling process to decide which generators will be
committed, may be affected as more Synchronous Generators (SGs) could be
dispatched for the sole purpose of providing inertia.

Several papers have approached the problem of including inertia requirements
in the UC problem. In [5, 6, 7], the authors use the swing equation of Center-
of-Inertia (CoI), which allows them to derive the RoCoF constraint and study
its effect on the UC schedule. However, this approach oversimplifies the
problem as it neglects metrics related to frequency deviation from the
setpoint. This problem was addressed in [8] and [9] with the inclusion of a
constraint limiting the post-disturbance maximum frequency deviation (i.e.,
frequency nadir). In [8], the analytic form of frequency nadir as a function
of active power disturbance is derived using a system frequency model obtained
from [10]. The nadir expression is then linearized and added to the UC model,
while considering a fixed sudden load increase. On the other hand, the authors
of [9] bypass the explicit modeling of turbine and governor control, as well
as their impact on frequency dynamics, by imposing strict assumptions on
system damping and total frequency response provision at each node. Moreover,
they look at the impact of wind uncertainty on inertia requirements. While the
simplifications proposed in [9] enable the formulation of a nadir constraint
without the explicit consideration of second-order frequency dynamics, they
oversimplify the actual control implementation and disregard the aggregate
impact of governor damping. Furthermore, none of the aforementioned studies
incorporate the converter interface of RES and the impact of respective
control schemes on the UC formulation.

This study builds on the work of [8] and [9], and improves on it in several
ways. First, we improve the frequency dynamics model in [8] by including the
state-of-the-art converter control schemes of inverter-based generation, more
specifically Virtual Synchronous Machine (VSM) and droop control. In contrast
to the existing literature, where SG inertia and damping constants are usually
numerically modified in order to compensate for high RES integration, we
analyze a realistic model of a low-inertia system comprising both SG and
converter dynamic models. This allows us to derive detailed analytic
expressions of relevant frequency metrics as functions of multiple system
variables (e.g., inertia, damping, aggregate droop gain, etc.) to be
determined by the UC model, as opposed to the approach in [9] where inertia
constant was the only decision variable of interest. Moreover, in addition to
frequency nadir and RoCoF, we incorporate the limitation on quasi steady-state
frequency deviation into the UC formulation. Secondly, a more straightforward
method is proposed to extract bounds for decision variables contributing to
frequency nadir, which allows us to incorporate the non-linear nadir
constraint in the UC problem in a more efficient way compared to [8].
Furthermore, similar to [9], this paper includes both the wind uncertainty and
potential loss of generation in the UC model. However, we present a more
comprehensive approach towards event probability computation and structuring
of the scenario tree for the two-stage stochastic UC problem.

The rest of the paper is structured as follows. In Section II, the derivation
of post-contingency frequency dynamics in a low-inertia, multi-machine system
is discussed. The obtained time-domain, analytic expressions are then
linearized and incorporated into a stochastic UC formulation in Section III.
Subsequently, the modeling of uncertainties, namely equipment failure and wind
power, in the form of probabilistic scenarios is presented in Section IV.
Section V provides the mathematical formulation of the stochastic UC problem.
Finally, Section VI presents and discusses the simulation results using a
modified version of the IEEE RTS-96 system, whereas Section VII draws the main
conclusions and discusses the outlook of the study.

##  II Low-Inertia System Frequency Dynamics

###  II-A Inertial Response and Primary Frequency Control Model

We first focus on deriving a simplified, but sufficiently accurate, uniform
frequency response model of a low-inertia system previously introduced in
[11]. Let us consider a system comprised of traditional $(i\in\mathcal{I})$
and converter-based $(j\in\mathcal{J})$ generators depicted in Fig. 1.

Figure 1: Uniform system frequency dynamics model.

The generator dynamics are described by the swing equation, with $M_{g}$ and
$D_{g}$ denoting the normalized inertia and damping constants corresponding to
the synchronous generators’ CoI. The low-order model proposed in [12] is used
for modeling the governor droop and turbine dynamics; $T_{i}$ are the turbine
time constants, $R_{i}$ and $K_{i}$ are the respective droop and mechanical
power gain factors, while $F_{i}$ refers to the fraction of total power
generated by the turbines of synchronous machines. Furthermore, we incorporate
the impact of grid-forming converters, as they are the only type of power
electronic-interfaced units providing frequency support [13, 14]. A particular
focus is set on droop $(d\in\mathcal{J}_{d}\subseteq\mathcal{J})$ and VSM
$(v\in\mathcal{J}_{v}\subseteq\mathcal{J})$ control schemes, as two of the
currently most prevalent emulation techniques in the literature, which in fact
have equivalent properties in the grid-forming mode of operation [15]. Here,
$T_{d}=T_{v}\equiv T_{j}$ are the time constants of all converters, $R_{d}$
and $K_{d}$ are the respective droop and electrical power gain factors,
whereas $M_{v}$ and $D_{v}$ denote the normalized virtual inertia and damping
constants of VSM converters.

###  II-B Analytic Derivation of Frequency Metrics

From Fig. 1 we can now derive a transfer function $G(s)$ of a general-order
system dynamics, as follows:

$\displaystyle G(s)$ $\displaystyle=\dfrac{\Delta f}{\Delta
P_{e}}=\bigg{(}\underbrace{(sM_{i}+D_{i})+\sum\limits_{i\in\mathcal{I}}\dfrac{K_{i}(1+sF_{i}T_{i})}{R_{i}(1+sT_{i})}}_{\text{traditional
generators}}$
$\displaystyle+\underbrace{\sum\limits_{d\in\mathcal{J}_{d}}\dfrac{K_{d}}{R_{d}(1+sT_{d})}}_{\text{droop
converters}}+\underbrace{\sum\limits_{v\in\mathcal{J}_{v}}\dfrac{sM_{v}+D_{v}}{1+sT_{v}}}_{\text{VSM
converters}}\bigg{)}^{-1}.$ (1)

Assuming similar time constants $(T_{i}\approx T)$ of all synchronous
machines, usually 2-3 orders of magnitude higher than the ones of converters,
justifies the approximation $T\gg T_{j}\approx 0$. Now we can transform (1)
into the following expression:

$G(s)=\frac{1}{MT}\frac{1+sT}{s^{2}+2\zeta\omega_{n}s+\omega_{n}^{2}},$ (2)

where the natural frequency $(\omega_{n})$ and damping ratio $(\zeta)$ are

$\omega_{n}=\sqrt{\frac{D+R_{g}}{MT}},\quad\zeta=\frac{M+T(D+F_{g})}{2\sqrt{MT(D+R_{g})}},$
(3)

and parameters $(M,D)$ and $(F_{g},R_{g})$ represent weighted system and
synchronous generator averages, respectively. More details on mathematical
formulation can be found in [11].

Assuming a stepwise disturbance in the electrical power $\Delta
P_{e}(s)=-\Delta P/s$, we can derive the time-domain expression for frequency
deviation as well as the nadir ($\dot{f}_{\mathrm{max}}$), RoCoF
($\dot{f}_{\mathrm{max}}$) and steady-state deviation ($\Delta
f_{\mathrm{ss}}$) frequency metrics:

$\displaystyle\Delta f_{\mathrm{max}}$ $\displaystyle=-\frac{\Delta
P}{D+R_{g}}\left(1+\sqrt{\dfrac{T(R_{g}-F_{g})}{M}}e^{-\zeta\omega_{n}t_{m}}\right),$
(4a) $\displaystyle\dot{f}_{\mathrm{max}}$
$\displaystyle=\dot{f}(t_{0}^{+})=-\frac{\Delta P}{M},$ (4b)
$\displaystyle\Delta f_{\mathrm{ss}}$ $\displaystyle=-\frac{\Delta
P}{D+R_{g}},$ (4c)

with the introduction of the new variable
$\omega_{d}=\omega_{n}\sqrt{1-\zeta^{2}}$.

The accuracy of the proposed model has already been investigated and verified
in [11]. We can conclude that the frequency metrics of interest are directly
dependent on the average system parameters $M$, $D$, $R_{g}$ and $F_{g}$, and
thus they could be regulated through the UC model. In particular, RoCoF and
steady-state deviation can be explicitly controlled via
$\dot{f}_{\mathrm{max}}\sim M^{-1}$ and $\Delta
f_{\mathrm{ss}}\sim(D+R_{g})^{-1}$, while nadir can be modeled using a highly
non-linear function $\Delta f_{\mathrm{max}}\left(M,D,R_{g},F_{g}\right)$.

##  III Formulation of Frequency Constraints

The aforementioned frequency expressions in (4) are incorporated as
constraints into the stochastic UC problem, converted into SI and bounded by
prescribed ENTSO-e thresholds [16], as follows:

$\displaystyle\left|\frac{f_{b}\Delta
P}{D+R_{g}}\left(1+\sqrt{\dfrac{T(R_{g}-F_{g})}{M}}e^{-\zeta\omega_{n}t_{m}}\right)\right|$
$\displaystyle\leq\Delta f_{\mathrm{lim}},$ (5a)
$\displaystyle\left|\frac{f_{b}\Delta P}{M}\right|\leq\dot{f}_{\mathrm{lim}},$
(5b) $\displaystyle\left|\frac{f_{b}\Delta P}{D+R_{g}}\right|\leq\Delta
f_{\mathrm{ss,lim}},$ (5c)

with $f_{b}=50\,\mathrm{Hz}$ being the base frequency; $\Delta
f_{\mathrm{lim}}=0.4\,\mathrm{Hz}$ is the Under-Frequency Load Shedding (UFLS)
trigger, while $\dot{f}_{\mathrm{lim}}=0.5\,\mathrm{Hz}/\mathrm{s}$ and
$\Delta f_{\mathrm{ss,lim}}=0.2\,\mathrm{Hz}$ are the maximum permissible
RoCoF and steady-state frequency deviation.

Constraints (5b) and (5c) are linear, unlike the non-linear frequency nadir
constraint (5a). In order to avoid the high computational burden of a Mixed-
Integer Non-Linear Program (MINLP) formulation and have a measurable
optimality gap, we use a linear approximation of (5a) which allows us to
maintain a Mixed-Ineger Linear Program (MILP) formulation of the stochastic UC
problem.

###  III-A Piece-wise Linearization of Nadir Expression

The study in [8] proposes a Piece-Wise Linearization (PWL) technique for
obtaining a linearized expression for frequency nadir in order to subsequently
integrate it into a UC problem. To improve clarity, this technique is outlined
here and in Appendix A before comparing its computational burden against our
proposed approach introduced in Section III-B. Let us recall from Section II-B
that the frequency nadir expression is a function of four variables
$(R_{g},F_{g},M,D)$, and as such too complicated to be directly handled by the
PWL. Considering that the aggregate damping constant is of the form
$D(D_{i},D_{v},R_{d})$, with respective damping and droop gains usually
strictly prescribed within narrow ranges by the system operator, it is
justifiable to assume a constant $D$ and therefore becomes $\Delta
f_{\mathrm{max}}(M,R_{g},F_{g})$. Hence, the PWL formulation aims to minimize
the following objective function

$\displaystyle\min_{\Psi}\sum_{\eta}\bigg{(}$
$\displaystyle\max_{1\leq\nu\leq\overline{\nu}}\left\\{a_{\nu}R_{g}^{(\eta)}+b_{\nu}F_{g}^{(\eta)}+c_{\nu}M^{(\eta)}+d_{\nu}\right\\}$
$\displaystyle-\Delta
f_{\mathrm{max}}\left(R_{g}^{(\eta)},F_{g}^{(\eta)},M^{(\eta)}\right)\bigg{)}^{2}$
(6)

with $\Psi=\\{a_{\nu},b_{\nu},c_{\nu},d_{\nu},\forall\nu\\}$, being the set of
optimization variables, $\eta$ denoting the evaluation point and $\nu$
referring to the number of PWL segments. The objective function (6) penalizes
the difference between the appropriate PWL segment and the nadir function at
all evaluation points. Given the convex nature of the nadir function, the
inner $\max$ operator chooses the appropriate PWL segment for each evaluation
point by looking at which segment is closest to the curve at that specific
point. To improve understanding, we provide an illustration of the PWL method
from [8] on a simple one variable function in Appendix A.

Upon obtaining the optimal solution of the model in (6), denoted as
$(a^{*}_{\nu},b^{*}_{\nu},c^{*}_{\nu},d^{*}_{\nu})$, the nadir constraint can
be integrated into the MILP UC model by adding a set of inequalities described
in Appendix A, along with the nadir threshold constraint of the form
$f_{b}\,t_{3}\leq\Delta f_{\mathrm{lim}}$. The results for the approximation
of frequency nadir function for a test system of 20 generators described in
Section VI are shown in Fig. 2, where a loss of the largest unit is
considered. Note that Fig. 2 showcases the surface plot for a fixed inertia
constant $M$ and thus ignores one degree of freedom. The original surface is
presented in blue, whereas its PWL-approximation segments are the planes
depicted in various colors. It is important to note that the optimization
problem (6) is computationally intensive and thus in order to obtain results
within reasonable computational time, the number of PWL segments used for the
approximation as well as the number of evaluation points have to be kept low.

Figure 2: PWL of the nadir constraint for $M=9$.

###  III-B Extracting Bounds on Relevant Variables

Figure 3: All possible values of the nadir after a generator loss.

An alternative approach for linearizing the nadir constraint and integrating
it into the UC problem is to confine the values of $R_{g}$, $F_{g}$, $M$ and
$D$ within a plausible range to guarantee that the nadir threshold in (5a) is
not violated. With this approach the damping variable $D$ can easily be
included and does not need to be set constant. The scatter plot presented in
Fig. 3 reflects all possible values of frequency nadir after the loss of the
largest generator, for the same system as in Section III-A. In the general
case, for a system that comprises $|\mathcal{I}|$ generators, there will be
$2^{|\mathcal{I}|-1}$ possible generator commitment combinations after a
generator outage. By obtaining the set of these dispatch combinations, the
values of $R_{g}$, $F_{g}$ and $M$ at which the UFLS threshold is not violated
can be extracted, corresponding to the points below the shaded plane in Fig.
3. Subsequently, these values are used to substitute the nadir constraint in
the unit commitment as follows:

$F_{g}\geq F_{g}^{\mathrm{lim}},\quad R_{g}\geq R_{g}^{\mathrm{lim}},\quad
M\geq M^{\mathrm{lim}}.$ (7)

Table I provides a comparison of the proposed method to the PWL technique, in
terms of computational time that is needed to obtain the equivalent linear
nadir equations for a single value of $\Delta P$. It is clear that the PWL is
more time intensive, especially when aiming for an increased precision. It
should be noted though that for very large systems, the calculations of
$2^{|\mathcal{I}|-1}$ combinations for the bound extraction method would
become more computationally expensive. Indeed, the computation time increases
by a factor of $2^{\Delta|\mathcal{I}|}$ for every additional
$\Delta|\mathcal{I}|$ generators included in the system. For the purposes of
this paper, the proposed bound extraction method will be used as it is
significantly faster and introduces less error when applied to the
20-generator test system under investigation.

TABLE I: Computational cost of the linearization methods.

Linearization technique | Computational time [s]
---|---
PWL ($m=3$, $k=4$) | 70
PWL ($m=4$, $k=4$) | 7200
Bound extraction | 20

##  IV Modeling equipment-failure and wind power uncertainties

This section describes the modeling of uncertainty pertaining to equipment
failure and wind power production during power system operation. The uncertain
nature of wind power production is modeled using a set of scenarios $\Omega$
that captures the spatio-temporal interdependence of forecast errors, for
every wind farm location and during the whole scheduling horizon. Each wind
power realization scenario $\omega$ has the same probability of occurrence
denoted as $\pi_{\omega}$.

In terms of equipment failure uncertainty, we consider as the set of credible
contingencies, the unforeseen outages of synchronous generators, whereas
transmission assets are assumed to be $100\,\%$ reliable. In order to reduce
the computational burden, we follow the assumption from [17] considering that
the generation outages happen at a discrete time period, while failed assets
remain unavailable for the rest of the scheduling horizon, i.e., the Mean Time
To Repair (MTTR) is greater than the scheduling horizon of the day-ahead
electricity market. For the purpose of assessing the impact of frequency
constraints on the unit commitment schedule, we consider as contingency period
the one in which the power system faces the highest wind power penetration as
scarcity of inertia is most likely to occur in this time due to the
displacement of synchronous generators from the day-ahead schedule.

To calculate the probability $\pi_{c}$ associated with contingency scenario
$c$, we index the set of credible contingencies by
$\kappa=1,2,\dots,\mathcal{K}$. We denote by $A(\kappa,\tau)$ the random event
of contingency $\kappa$ happening within time period $\tau$. Random event
$B(\kappa)$ corresponds to contingency $\kappa$ not occurring during the
entire scheduling horizon. We further denote as $\lambda_{\kappa}$ the inverse
of the Mean Time To Failure $(\mathrm{MTTF})$ $\kappa$, i.e.,
$\lambda_{\kappa}=1/\mathrm{MTTF}$. Considering we are looking at only one
hour in which the outages may occur, the probability $\pi_{c}$ for each
contingency scenario $\kappa$ is derived from the probabilities of occurrence
of random events $A(\kappa,\tau)$ and $B(\kappa)$ that are calculated using
the following expressions according to [18]:

$\displaystyle\pi[A(1,\tau)]=\mathrm{exp}(-\lambda_{1}\tau)(\mathrm{exp}(\lambda_{1})\\!-\\!1),$
(8a)
$\displaystyle\pi[A(1,\tau)]=\pi[A(2,\tau)]=\ldots=\pi[A(\mathcal{K},\tau)],$
(8b) $\displaystyle\pi[B(1)]=\mathrm{exp}(-\lambda_{1}\tau),$ (8c)
$\displaystyle\pi[B(1)]=\pi[B(2)]=\ldots=\pi[B(\mathcal{K})].$ (8d)

Assuming statistical independence between all contingencies, the probability
$\pi_{c_{0}}$ of the no-contingency scenario is equal to

$\displaystyle\pi_{c_{0}}$
$\displaystyle=\prod_{\kappa=1}^{\mathcal{K}}\pi[B(\kappa)],$ (9)

while the probability $\pi_{c_{\kappa}}$ of losing a generator is equal to

$\displaystyle\pi_{c_{\kappa}}=\pi[A(1,\tau)]\prod_{\begin{subarray}{c}y=1\\\
y\neq
k\end{subarray}}^{\mathcal{K}}\pi[B(y)],\quad\forall\kappa=1,...,\mathcal{K}.$
(10)

It should be noted that the sum of probabilities $\pi_{c_{0}}$ and
$\pi_{c_{\kappa}}$ is lower than 1, since sequential contingencies are not
considered. For instance, setting $\mathrm{MTTF}$ equal to $1000\,\mathrm{h}$
for all generators, we obtain from (8)
$\pi[A(\kappa,\tau)]=0.9995\\!\times\\!10^{-3}$ and $\pi[B(\kappa)]=0.9990$.
According to (9) and (10) we obtain respectively $\pi_{c_{0}}=0.9960$ and
$\pi_{c_{\kappa}}=0.9965\times 10^{-3}$ and thus
$\sum_{\kappa=0}^{4}\pi_{c_{\kappa}}=0.9999\approx 1$.

Figure 4: Scenario tree for the two-stage stochastic UC problem.

Combining the scenarios modeling the equipment failure and wind power
uncertainty into a single scenario set $\mathcal{E}$, we define each scenario
$\xi$ as a pair of contingency $c$ and wind power realization $\omega$. For
each $\xi=\\{c,\omega\\}$ the corresponding probability of occurrence is given
as $\pi_{\xi}=\pi_{\omega}\cdot\pi_{c}$ and
$\sum_{\xi\in\mathcal{E}}\pi_{\xi}\approx 1$, assuming that equipment outages
and wind power production are statistically independent events. The structure
of the scenario set $\mathcal{E}$ used in the stochastic UC formulation is
illustrated as the scenario tree shown in Fig. 4, for $|\mathcal{K}|$
contingencies and $|\mathcal{W}|$ wind power scenarios.

##  V Stochastic Unit Commitment

This section provides the mathematical formulation of the stochastic unit
commitment [19], with an addition of frequency-related constraints. The
proposed model is a two-stage stochastic optimization problem which can be
written as:

$\displaystyle\min_{\Phi}\,\,\sum_{t\in\mathcal{T}}\sum_{i\in\mathcal{I}}\big{(}C^{\mathrm{SU}}_{i}y_{it}+C^{\mathrm{SD}}_{i}z_{it}+C_{i}p_{it}\big{)}\
+$ (11a)
$\displaystyle\sum_{t\in\mathcal{T}}\sum_{\xi\in\mathcal{E}}\pi_{\xi}\Big{[}\sum_{i\in\mathcal{I}}\big{(}C^{+}_{i}r^{+}_{i\xi
t}-C^{-}_{i}r^{-}_{i\xi
t}\big{)}+\sum_{n\in\mathcal{N}}C^{\mathrm{sh}}\,l^{\mathrm{shed}}_{n\xi
t}\Big{]}$ subject to
$\displaystyle\sum_{i\in\mathcal{I}_{n}}p_{it}\\!+\\!\sum_{j\in\mathcal{J}_{n}}w_{jt}-D_{nt}-$
$\displaystyle\quad\sum_{m:(n,m)\in\mathcal{L}}B_{nm}(\hat{\delta}_{nt}-\hat{\delta}_{mt})=0,\quad\forall
n,t,$ (11b) $\displaystyle
B_{nm}(\hat{\delta}_{nt}-\hat{\delta}_{mt})\leq\overline{f}_{nm},\quad\forall(n,m)\in\mathcal{L},t,$
(11c) $\displaystyle u_{it}-u_{i(t-1)}\leq u_{i\tau_{i}^{1}},\quad\forall
i,t,$ (11d) $\displaystyle u_{i(t-1)}-u_{it}\leq
1-u_{i\tau_{i}^{0}},\quad\forall i,t,$ (11e) $\displaystyle y_{it}\geq
u_{it}-u_{i(t-1)},\quad\forall i,t,$ (11f) $\displaystyle z_{it}\geq
u_{i(t-1)}-u_{it},\quad\forall i,t,$ (11g)
$\displaystyle\sum_{i\in\mathcal{I}_{n}}\Big{[}r^{+}_{i\xi t}-r^{-}_{i\xi
t}-p_{it}(1-\alpha_{i\xi t})\Big{]}$
$\displaystyle+\sum_{m:(n,m)\in\mathcal{L}}B_{nm}(\hat{\delta}_{nt}-\tilde{\delta}_{n\xi
t}-\hat{\delta}_{mt}+\tilde{\delta}_{m\xi t})$ (11h)
$\displaystyle+\sum_{j\in\mathcal{J}_{n}}(W^{*}_{j\xi
t}-w_{jt}-w^{\mathrm{spill}}_{j\xi t})+l^{\mathrm{shed}}_{n\xi
t}=0,\quad\forall n,\xi,t,$ $\displaystyle p_{it}+r^{+}_{i\xi
t}\leq\overline{P_{i}}u_{it},\quad\forall i,\xi,t,$ (11i) $\displaystyle
p_{it}-r^{-}_{i\xi t}\geq\underline{P_{i}}u_{it},\quad\forall i,\xi,t,$ (11j)
$\displaystyle p_{it}-p_{i(t-1)}+r^{+}_{i\xi t}-r^{+}_{i\xi(t-1)}\leq
R^{\mathrm{U}}_{i},\quad\forall i,\xi,t,$ (11k) $\displaystyle
p_{it}\\!-\\!p_{i(t-1)}-r^{-}_{i\xi
t}\\!+\\!r^{-}_{i\xi(t-1)}\geq-R^{\mathrm{D}}_{i},\quad\forall i,\xi,t,$ (11l)
$\displaystyle r^{+}_{i\xi t}\leq R_{i}^{+}\alpha_{i\xi t},\quad\forall
i,\xi,t,$ (11m) $\displaystyle r^{-}_{i\xi t}\leq R_{i}^{-}\alpha_{i\xi
t},\quad\forall i,\xi,t,$ (11n) $\displaystyle
B_{nm}(\tilde{\delta}_{nt}-\tilde{\delta}_{mt})\leq\overline{f}_{nm},\quad\forall(n,m)\in\Lambda,t,$
(11o) $\displaystyle w^{\mathrm{spill}}_{j\xi t}\leq w_{j\xi t},\quad\forall
j,\xi,t,$ (11p) $\displaystyle l^{\mathrm{shed}}_{n\xi t}\leq
D_{nt},\quad\forall n,\xi,t,$ (11q) $\displaystyle k_{i\xi
t}=\frac{P_{i}K_{i}}{\sum_{i\in\mathcal{I}}P_{i}+W_{v}+W_{d}}u_{it}\alpha_{i\xi
t},\quad\forall i,\xi,t,$ (11r) $\displaystyle F_{\xi
t}=\sum_{i\in\mathcal{I}}\frac{F_{i}k_{i\xi t}}{R_{i}},\quad\forall t,\,\xi,$
(11s) $\displaystyle R_{\xi t}=\sum_{i\in\mathcal{I}}\frac{k_{i\xi
t}}{R_{i}},\quad\forall t,\,\xi,$ (11t) $\displaystyle M_{\xi
t}=\sum_{i\in\mathcal{I}}2H_{g}k_{i\xi t},\quad\forall t,\,\xi,$ (11u)
$\displaystyle\frac{\dot{f}_{\mathrm{lim}}}{f_{b}}(M_{\xi t}+M_{v})\geq\Delta
P_{\xi t},\quad\forall t,\,\xi,$ (11v) $\displaystyle F_{\xi t}\geq F_{\xi
t}^{\mathrm{lim}},\ R_{\xi t}\geq R_{\xi t}^{\mathrm{lim}},\ M_{\xi
t}+M_{v}\geq M_{\xi t}^{\mathrm{lim}},\quad\forall t,\,\xi,$ (11w)
$\displaystyle\frac{\Delta f_{\mathrm{ss,lim}}}{f_{b}}(D+R_{\xi t})\geq\Delta
P_{\xi t},\quad\forall t,\,\xi,$ (11x) $\displaystyle p_{it}\geq 0,\forall
i,t;\,w_{jt}\geq 0,\forall j,t;\,\hat{\delta}_{nt}\geq 0,\forall n,t;\,k_{i\xi
t}\geq 0,\forall i,\xi,t;$ $\displaystyle r^{+}_{i\xi t},r^{-}_{i\xi t}\geq
0,\forall i,\xi,t;\,l^{\mathrm{shed}}_{n\xi t}\geq 0,\forall
n,\xi,t;\,w^{\mathrm{spill}}_{j\xi t}\geq 0,\;\forall j,\xi,t;$ $\displaystyle
F_{\xi t},R_{\xi t},M_{\xi t}\geq
0,\forall\xi,t;\,u_{it},y_{it},z_{it}\in\\{0,1\\},$ (11y)

where $\Phi=\\{p_{it},\,u_{it},\,y_{it},\,z_{it},\,\forall
i,t;\;w_{jt},\forall j,t;\;\hat{\delta}_{nt},\forall n,t;$
$\tilde{\delta}_{n\xi t},\forall n,\xi,t;\,r^{+}_{i\xi t},\,r^{-}_{i\xi
t},\forall i,\xi,t;\,w^{\mathrm{spill}}_{j\xi t},\forall
j,\xi,t;\,l^{\mathrm{shed}}_{n\xi t},\forall n,\xi,t;\,$ $k_{i\xi t},\forall
i,\xi,t;\,F_{\xi t},\,R_{\xi t},\,M_{\xi t},\forall\xi,t\\}$ is the set of
optimization variables.

The objective function (11a) to be minimized is the total expected system cost
that comprises the day-ahead energy and the real-time balancing costs. The
day-ahead component consists of the fuel costs $C_{i}$ as well as the start-up
and shut-down costs. The real-time component includes the re-dispatch cost
from the deployment of upward and downward reserves based on the corresponding
offer prices $C^{+}_{i}$ and $C^{-}_{i}$, as well as the involuntary load
shedding at the value of lost load $C^{\mathrm{sh}}$.

Equation (11b) enforces the nodal power balance of the day-ahead schedule,
while network power flows at the day-ahead stage are restricted by the
transmission capacity limits in (11c). Constraints (11d)-(11e) model the
minimum online and offline time of conventional units based on commitment
variable $u_{it}$, where parameters $\tau^{1}_{i}$ and $\tau^{0}_{i}$ are
defined as $\tau_{i}^{1}=\min\\{t+T^{1}_{i}-1,\ T\\}$ and
$\tau_{i}^{0}=\min\\{t+T^{0}_{i}-1,\ T\\}$, and $T^{1}_{i}$ and $T^{0}_{i}$
denote the duration that unit $i$ should remain online and offline,
respectively. Constraints (11f)-(11g) model the start-up and shut-down of
conventional units using the binary variables $y_{it}$ and $z_{it}$,
respectively. The real-time power balance for every uncertainty realization
$\xi$ is enforced by constraint (11h). Parameter $\alpha_{i\xi t}$ models the
availability of the generators to provide reserves, i.e., $\alpha_{i\xi t}$ is
equal to 1 if generator $i$ at scenario $\xi$ and time $t$ is online and able
to provide reserves and zero otherwise. The scheduled energy production and
the deployment of upward ($r^{+}_{i\xi t}$) and downward ($r^{-}_{i\xi t}$)
reserves in each scenario $\xi$ are bounded by the generation capacity limits
of each unit by constraints (11i)-(11j), whereas constraints (11k)-(11l)
enforce the upward and downward ramping limits accounting for the real-time
reserve activation. Constraints (11m)-(11n) account for the limits of reserve
capacity offers. Transmission capacity limits during real-time operation are
enforced by constraint (11o), whereas wind spillage $w^{\mathrm{spill}}_{j\xi
t}$ and load shedding $l^{\mathrm{shed}}_{n\xi t}$ are bounded by the wind
power realization and the nodal demand through constraints (11p) and (11q),
respectively.

The set of constraints (11r)-(11w) models the frequency limits of the power
system. The equality constraint (11r) defines $k_{i\xi t}$ as the gain factor
$K_{i}$ of generator $i$ scaled by the ratio of its capacity over the total
system capacity, which in turn is multiplied by the binary variable $u_{it}$
and the parameter $\alpha_{i\xi t}$ to indicate that a unit can only provide
inertial response if it is committed and does not face an outage. Similarly,
constraints (11s)-(11u) define average system variables for power fraction,
droop and inertia, respectively. Constraint (11v) enforces the RoCoF limit,
while nadir equivalent and quasi steady-state frequency bounds are imposed by
constraints (11w) and (11x). Finally, constraints (11y) are variable
declarations.

##  VI Case Study

###  VI-A System Description

In order to analyze the performance of the stochastic UC model presented in
Section V, we investigate a modified version of the IEEE RTS-96 power system
from [20] depicted in Fig. 5, with 48 buses comprising areas 1 and 2 of the
original system. Table II shows the relevant parameters of different thermal
plant types. The studied system includes 20 generators and 16 wind farms. We
assume that six wind farms are providing virtual inertia; four of them via VSM
control, and the remaining two through equivalent droop regulation. The UC is
ran for two days without frequency constraints in order to initialize the
system prior to introducing the frequency constraints on days 3, 4 and 5. This
is done to ensure the impact of start-up costs are well distributed and not
concentrated on one day. Therefore, the total simulation horizon is five days
($T=120\,\mathrm{h}$) whereas the UC schedule is optimized separately for each
day, with the last hour of each day used as an input for the next.

Figure 5: Modified IEEE RTS-96 system diagram comprised of 2 areas, 16 wind
farms and 20 synchronous generators [20].

As the set of possible contingencies we consider the failure of synchronous
generators $i=\\{1,6,8,10\\}$. These generators are of various capacities,
ranging from the smallest to the largest unit in the system. The hour 19 of
day 3 (i.e., $t=67\,\mathrm{h}$) is selected to be the time instance of a
possible contingency, as this is the hour with high wind penetration and low
demand. Ten wind power scenarios are considered, which brings the total number
of scenarios to 50. The optimization problem is formulated in Python and uses
the Gurobi solver with default parameterization.

TABLE II: Parameters of the thermal plants and VSM.

Type | $H_{g}\,[\mathrm{s}]$ | $K_{g}\,[\mathrm{p.u.}]$ | $F_{g}\,[\mathrm{p.u.}]$ | $R_{g}\,[\mathrm{p.u.}]$ | $D_{g}\,[\mathrm{p.u.}]$
---|---|---|---|---|---
Nuclear | 4.5 | 0.98 | 0.25 | 0.04 | 0.6
CCGT | 7.0 | 1.1 | 0.15 | 0.01 | 0.6
OCGT | 5.5 | 0.95 | 0.35 | 0.03 | 0.6
VSM | 6.0 | 1.0 | - | - | 0.6
Droop | - | 1.0 | - | 0.05 | -

###  VI-B Results

In this section, the simulation results from the stochastic UC are presented.
Fig. 6 showcases the load and wind power profiles as well as the aggregate
dispatch of synchronous generation for two UC runs: (i) without frequency
constraints; and (ii) with frequency constraints. Furthermore, Table III
indicates the difference in the total number of generators committed between
the two runs. Both Fig. 6 and Table III suggest that, although the amount of
committed generators increases significantly, the total SG production is only
slightly changed. This is justified by the fact that the extra generators are
solely committed for the purpose of providing inertia, and are thus operating
at their technical minimum. The production surplus arising from the
additionally committed units is compensated by wind curtailment and other
generators reducing their power output.

Figure 6: UC dispatch of synchronous generation for respective load and wind
profiles. TABLE III: Comparison of the total number of dispatched generators
through UC for each hour.

Hour | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73
---|---|---|---|---|---|---|---|---|---
w/o FC | 6 | 5 | 4 | 4 | 4 | 4 | 4 | 4 | 4
w/ FC | 6 | 5 | 10 | 10 | 10 | 10 | 10 | 10 | 10

Figure 7: Impact of frequency constraints on the aggregate level of system
inertia.

Moreover, the evolution of aggregate system inertia over the course of the
whole scheduling horizon is depicted in Fig. 7. A noticeable step change in
total system inertia at hour $67$ reflects the violation of frequency
constraints under contingency, which subsequently triggers a dispatch of
auxiliary synchronous generators. While the inertia levels do not
differentiate between the two scenarios during the first two days, on the days
following the potential outage some carryover impacts can be observed. This is
a consequence of the commitment schedule being radically changed at hour
$t=67$, thus affecting the UC schedules in the following days.

Some insightful conclusions can be drawn from Fig. 8, where the difference
between the actual values of the frequency metrics and the respective ENTSO-E
thresholds are depicted. For this purpose, we define a constraint gap $\eta$
as a measure of the relative constraint distance to its limit, e.g.,
$\eta_{\mathrm{nadir}}=\Delta f_{\mathrm{max}}/\Delta f_{\mathrm{lim}}-1$.
After the completion of unit commitment, the constraints are re-evaluated
using the obtained $F_{g}$, $R_{g}$ and $M$ values in order to determine which
frequency criteria becomes binding at the instance of the fault. A negative
constraint gap corresponds to the non-binding constraint, i.e., the specific
frequency criteria being met. It should be noted that the positive values of
$\eta$ for $t<67\,\mathrm{h}$ indicate that the frequency threshold would be
violated if a fault occurs. No action is required however, considering that in
this case study we assume that the contingency can only occur in hour $t=67$.
The Fig. 8 indicates that without explicit frequency metric constraints all of
these constraints would be violated. Moreover, it can be observed that when
including the frequency constraints, the RoCoF constraint is closest to its
limit - corresponding to the smallest constraint gap - and thus binding. The
constraint gap difference between the two scenarios at hour $t=67$ clearly
highlights the importance of including the frequency constraints in UC in
order to avoid large frequency excursions and undesired triggering of
protection and UFLS schemes. The same observations are also reflected in Fig.
9 through time-domain frequency response of the system. Understandably, the
values of RoCoF, nadir and steady-state deviation are reduced compared to the
scenario without frequency constraints, such that all of the ENTSO-E criteria
are fulfilled.

Figure 8: Constraint gaps for different frequency metrics. Dashed lines refer
to the scenario without frequency constraints.

Figure 9: Frequency evolution of the CoI with and without FC for 20 seconds
after the fault instance at $t=67\,\mathrm{h}$.

Finally, we investigate the economic impact of including the frequency
constraints into the stochastic UC model. The breakdown of operational costs
for day 3 is presented in Table IV. The addition of frequency constraints
leads to a $5\,\%$ increase in total expected system costs and a significant
increase in start-up costs by $185\,\%$. This is due to six extra generators
being turned on for providing inertia at the period of a potential generation
failure, as shown in Table III. A large increase is also seen in reserve
scheduling costs, as the reserves are now not only scheduled to cover wind
power uncertainty but also for possible contingencies. The change in cost is
of course highly dependent on the specific system and the considered
contingencies.

##  VII Conclusion

This paper includes frequency constraints in the UC problem of a system with
large wind power penetration in order to investigate the impact of frequency
dynamics on unit scheduling. By employing the analytic expressions for post-
contingency frequency response of a multi-machine system, we define a set of
constraints reflecting the frequency nadir, RoCoF and quasi steady-state
deviation. The highly non-linear frequency nadir constraint is linearized
using two approaches: (i) a PWL technique adapted from the literature; and
(ii) a proposed simple and efficient method for extracting bounds on decision
variables of interest, which is shown to be computationally superior to PWL.
Using the latter approach, the stochastic UC problem is formulated as an MILP,
with an objective of minimizing the expected system costs against wind power
production and generation outage uncertainties.

TABLE IV: Unit commitment costs $\mathrm{[}\$]$ breakdown on day 3.

Scenario | Total costs | Start-up | Operation | Reserves
---|---|---|---|---
w/o FC | $410\,545$ | $788$ | $405\,561$ | $4\,196$
w/ FC | $432\,383$ | $2\,248$ | $420\,886$ | $9\,249$
Difference | $5.32\,\%$ | $185\,\%$ | $3.78\,\%$ | $120\,\%$

Our results show that the inclusion of frequency constraints in the UC model
significantly affects the dispatch of synchronous generators and consequently
the expected system costs. Indeed, during anticipated critical events such as
the loss of generation, additional synchronous machines are needed for
providing sufficient inertia and damping in the process of frequency
containment. Such actions lead to a drastic increase in the UC costs,
especially start-up and reserve scheduling, which poses a new challenge as the
operator must find a way to remunerate the units committed for the sole
purpose of frequency regulation. This is an exciting avenue for future work.

## Appendix A

For clarity, we visually illustrate in Fig. 10 how the PWL optimization
problem from (6) is solved. In this example there are four evaluation points
set at -7.5, -2.5, 2.5 and 7.5, and the respective function is approximated
with four segments. At each evaluation point the model identifies the segment
closest to the original curve, and subsequently aims to minimize the overall
shaded area. The proposed technique can then be expanded and employed on a
function of three variables, as we have done with the expression for frequency
nadir.

Figure 10: Illustration of the PWL method on a 2-D function.

References section

## References

  * [1] F. Milano, F. Dörfler, G. Hug, D. J. Hill, and G. Verbič, “Foundations and challenges of low-inertia systems,” in _Power Systems Computation Conference (PSCC)_ , Jun 2018.
  * [2] P. Tielens and D. Van Hertem, “Grid inertia and frequency control in power systems with high penetration of renewables,” _Young Researchers Symposium in Electrical Power Engineering_ , 2012.
  * [3] J. Matevosyan and P. Du, “Wind integration in ERCOT,” in _Integration of Large-Scale Renewable Energy into Bulk Power Systems_. Springer, 2017, pp. 1–25.
  * [4] EirGrid and SONI, “DS3: System services review TSO recommendations,” EirGrid, Tech. Rep., May 2013.
  * [5] E. Davarinejad, M. R. Hesamzadeh, and H. Chavez, “Incorporating inertia constraints into the power market,” Energiforsk, Tech. Rep., 2017.
  * [6] F. Pérez-Illanes, E. Álvarez-Miranda, C. Rahmann, and C. Campos-Valdés, “Robust unit commitment including frequency stability constraints,” _Energies_ , vol. 9, no. 11, p. 957, 2016.
  * [7] P. Daly, D. Flynn, and N. Cunniffe, “Inertia considerations within unit commitment and economic dispatch for systems with high non-synchronous penetrations,” in _PowerTech, 2015 IEEE Eindhoven_ , 2015.
  * [8] H. Ahmadi and H. Ghasemi, “Security-constrained unit commitment with linearized system frequency limit constraints,” _IEEE Trans. Power Syst._ , vol. 29, pp. 1–11, 07 2014.
  * [9] F. Teng, V. Trovato, and G. Strbac, “Stochastic scheduling with inertia-dependent fast frequency response requirements,” _2016 IEEE Power and Energy Society General Meeting (PESGM)_ , 2016.
  * [10] D. Lee Hau Aik, “A general-order system frequency response model incorporating load shedding: Analytic modeling and applications,” _IEEE Trans. Power Syst._ , vol. 21, pp. 709 – 717, 06 2006.
  * [11] U. Markovic, Z. Chu, P.Aristidou, and G. Hug, “LQR-based adaptive virtual synchronous machine for power systems with high inverter penetration,” _IEEE Trans. Sustainable Energy_ , Jan 2019 (forthcoming).
  * [12] P. M. Anderson and M. Mirheydar, “A low-order system frequency response model,” _IEEE Trans. Power Syst._ , vol. 5, no. 3, pp. 720–729, Aug 1990\.
  * [13] J. Rocabert, A. Luna, F. Blaabjerg, and P. Rodríguez, “Control of power converters in AC microgrids,” _IEEE Trans. Power Electron._ , vol. 27, no. 11, pp. 4734–4749, Nov 2012.
  * [14] U. Markovic, O. Stanojev, P. Aristidou, and G. Hug, “Partial grid forming concept for 100% inverter-based transmission systems,” in _IEEE PES General Meeting_ , Aug 2018.
  * [15] U. Tamrakar, D. Shrestha, M. Manisha, B. P. Bhattarai, T. M. Hansen, and R. Tonkoski, _Virtual Inertia: Current Trends and Future Directions_. MDPI AG, June 2017.
  * [16] ENTSO-E, “Frequency stability evaluation criteria for the synchronous zone of continental europe,” ENTSO-E, Tech. Rep., 2016.
  * [17] A. Ahmadi-Khatir, A. J. Conejo, and R. Cherkaoui, “Multi-area energy and reserve dispatch under wind uncertainty and equipment failures,” _IEEE Trans. Power Syst._ , vol. 28, pp. 4373–4383, 2013.
  * [18] A. J. Conejo, M. Carrión, and J. M. Morales, _Decision making under uncertainty in electricity markets_. Springer, 2010, vol. 1.
  * [19] S. Delikaraoglou and P. Pinson, “High-quality wind power scenario forecasts for decision-making under uncertainty in power systems,” in _13th IEEE International Workshop on Large-Scale Integration of Wind Power and Transmission Networks (WIW 2014)_ , 2014.
  * [20] H. Pandzic, Y. Dvorkin, T. Qiu, Y. Wang, and D. Kirschen, Unit Commitment under Uncertainty - GAMS Models, Library of the Renewable Energy Analysis Lab (REAL), University of Washington, Seattle, USA. [Online]. Available at: http://www.ee.washington.edu/research/real/gams_code.html.


# Sensitivity of Photovoltaic Cells Efficiency to Initial Conditions in
Various Aggregation Designs

Baharak Mohamad Jafari Navadel, Esfandyar Faizi, Baharam Ahansaz ,
Jaber Jahanbin Sardroodi,
Physics Department, Azarbaijan Shahid Madani University, Tabriz, Iran
<EMAIL_ADDRESS>
###### Abstract

It is thought that nature already exploits quantum mechanical properties to
increase the efficiency of solar energy harvesting devices. So, the operation
of these devices can be enhanced by clever design of a nanoscopic, quantum
mechanical system where the quantum coherence plays a crucial role in this
process. In this investigation, we develop a donor-acceptor two level trap
dipole model converging the key role of quantum coherence and aggregation
effects along with different initial states. Our analysis reveals that
quenching unwanted emissions is achievable by preparing the system in specific
initial state under the effect of optimal spatial aggregation. Interestingly
it is observed that characterizing aggregation-induced properties and quantum
effects of bandgap engineering can increase the power enhancement up to
$35.87\%$ compared with classical counterparts. This encouraging trend
suggests a promising novel design aspect of nature-mimicking photovoltaic
devices.

Keywords: Photovoltaic cell, Quantum heat engine, Quantum coherence,
H-aggregate, J-aggregate, Dark states

##  I. INTRODUCTION

Photovoltaic cells (PV cells), also known as solar cells and photosynthesis,
akin to classical heat engines, convert solar energy directly into electrical
and chemical energy by the means of photovoltaic effect and a biological
process, respectively. The pursuit of enhancing the energy conversion
efficiency of PVs is the subject of extensive ongoing research. The
performance of a solar energy harvesting device can be improved through the
strategic design of a nanoscopic, quantum mechanical system. While
thermodynamic principles establish the well-known Shockley-Queisser efficiency
limit for classical photovoltaic devices, this limit can be surpassed by
deliberately utilizing quantum interference to disrupt the detailed balance
constraint [1]. It is hypothesized that nature leverages quantum mechanical
properties to enhance the light-harvesting efficiency of photosynthesis.
Prolonged quantum coherence has been observed in photosynthesis following
laser excitation [2, 3, 4, 5, 6, 7]. This observation has attracted
substantial interest in understanding how quantum coherence can be amplified
in complex biological environments and its potential crucial role in efficient
exciton transport processes [8, 9, 10, 11, 12]. The FMO complex which connects
the antenna to the reaction centre in the light harvesting apparatus of green
sulfur bacteria, stands as the most extensively investigated system within
this context [13]. Emulating photosynthesis offers a compelling avenue for
enhancing the efficiency of contemporary solar cell technology [14].
Establishing a connection between efficiency, functionality, and fortifying
room-temperature quantum effects in these nanoscale systems could profoundly
influence the design of future quantum-based nanotechnologies.

In a recent study, Dorfman and colleagues proposed an innovative solution to
enhance light reactions in photocells [15]. They analyzed these reactions as
quantum heat engines (QHEs). By treating the light-to-charge conversion as a
continuous Carnot-like cycle, they discovered that quantum coherence could
significantly increase the photocurrent in a photocell based on photosynthetic
reaction centers. Specifically, this boost amounted to at least $27\%$ when
compared to an equivalent classical photocell. In their theoretical framework,
the driving force behind this enhancement stems from the phenomenon of Fano
interference [16]. Fano interference has been experimentally shown to allow
optical systems to deviate from the thermodynamic detailed balance, which
typically constrains the efficiency of light-harvesting devices. Notably, this
violation of detailed balance was originally highlighted by Shockley and
Queisser in their seminal work on the fundamental limits of semiconductor
solar cells. In this regard, Scully and his colleagues theoretically
demonstrated that quantum coherence can enhance the performance of both a
solar cell and a photosynthetic reaction center [16, 17, 18]. Succeeding on
the work of Scully et al., Creatore et al. proposed a biologically inspired
photocell model enhanced by a delocalized dark quantum state involving two
dipole-coupled donors [19]. Next, the authors of Ref. [20] conducted an
investigation on the scenario involving three coupled donors. Furthermore,
Fruchtman et al. demonstrated that a photocell incorporating an asymmetric
pair of coupled chromophores can surpass the performance of those containing a
symmetric dimer or a pair of independent molecules [21]. The distinction
between coherent and incoherent energy transfer has been extensively
investigated in molecular crystals and aggregates. It is well-established that
the interaction between exciton coupling and energetic disorder regulates the
degree of exciton delocalization, which subsequently dictates the nature of
transport. As excitons become more delocalized, coherent effects become
increasingly significant.

In this paper we set out dipole coupled emitters and spatial nano-structure
two-level trap required to fine-tune the emission properties such as power
enhancement through preparing the system in different initial state and
specific aggregation model. It is well established that emission properties
are the reflection of interaction between optical molecular transition dipoles
in collective excitations. We host dipole molecules of solar materials in H
and J-aggregates producing delocalized states in order to supress
recombination and increase photo-generated carriers. Theoretical studies have
shown that H-aggregate is more advantageous for charge transport due to the
large $\pi-\pi$ overlap. Based on this scientific background, this research
has provided systematic insights into the correlation between aggregation
structure of dipole electron donors and emission induced properties causing.
More interestingly, the system preparation in multiple initial states made a
big improvement in power output enhancement intoducing unique, novel and more
enhanced photocell model.

##  II. MODEL DESCRIPTION

In this paper, we develop the concept of a photosynthesis-inspired paradigm as
QHE in order to increase the light-harvesting efficiency of quantum
photocells. The cyclic engine model considered here is a structure with three
effective sites which mimics the photosynthetic reaction centers apparatus as
in the primary model proposed by Dorfman [15]. According to crystallography
molecules in aggregates are most likely to be aligned collectively, not
independently boosting individual optoelectronic properties in tunable light-
matter interactions [22, 23, 24, 25, 26, 27, 28]. The well-defined aggregation
structure we consider here consists of two identical and initially uncoupled
donor chromophores, ($D$, e.g., polymer-based material) which flank an
acceptor molecule ($A$, e.g., fullerene-basedmaterial), as depicted in Figure
1. The optical excitation of donor configurations produced by solar radiation
is modeled as two-level systems with the ground state $|b\rangle$ and the
excited states $|a_{1}\rangle$ and $|a_{2}\rangle$. The exciton dynamics of
the donor aggregate structure is described by the Hamiltonian given by

$\displaystyle\mathcal{H}_{D}=\sum_{i=1,2}\hbar\omega_{i}\sigma_{i}^{+}\sigma_{i}^{-}+\sum_{i\neq
j}J_{ij}(\sigma_{i}^{+}\sigma_{j}^{-}+\sigma_{i}^{-}\sigma_{j}^{+}),$ (1)

where $\sigma_{i}^{+}=(\sigma_{i}^{-})^{\dagger}=|a_{i}\rangle\langle b|$ with
$i=1,2$ are the raising (lowering) operators. By considering the donor
molecules as a dipole moment with the optical transition dipole moment
$\mathbf{\mu_{i}}=e\langle a_{i}|\mathbf{r}|b\rangle$ $(i=1,2)$, the
electrostatic dipole-dipole coupling describes intermolecular interaction is
given by

$\displaystyle
J_{12}=\dfrac{1}{4\pi\varepsilon\varepsilon_{0}}\bigg{(}\frac{\mathbf{\mu_{1}}.\mathbf{\mu_{2}}}{|\mathbf{R_{12}}|^{3}}-3\frac{(\mathbf{\mu_{1}}.\mathbf{R_{12}})(\mathbf{\mu_{2}}.\mathbf{R_{12}})}{|\mathbf{R_{12}}|^{5}}\bigg{)}.$
(2)

Consider dipole moment $\mathbf{\mu_{1}}$ located at position $\mathbf{R_{1}}$
and dipole moment $\mathbf{\mu_{2}}$ at position $\mathbf{R_{2}}$ coupled to
electromagnetic radiation field through dipole transitions and the radius
vector from $\mathbf{\mu_{1}}$ to $\mathbf{\mu_{2}}$ is
$\mathbf{R_{12}}=\mathbf{R_{2}}-\mathbf{R_{1}}$. Typically, the strength of
the interaction term $J_{12}$ is much weaker than the optical excitation
energy and depends on the alignment of dipole moments. In a dimer system
studied by Creatore et al. [19], the donor dipole moment is always
perpendicular to the radius vector $\mathbf{R_{12}}$, resulting in the
vanishing of second term in Eq.(2). Thus, the electrostatic dipole-dipole
coupling is given by
$J_{12}\propto\mathbf{\mu_{1}}.\mathbf{\mu_{2}}=\mu_{1}\mu_{2}\mathrm{cos}\varphi$
where $\varphi$ is the angle between the two dipole moments. But with the
assumption that the donor dipole moments are parallel ($\varphi=0$) and when
the angle between them and the vertical axis denoted by $\theta$, the reduced
form of Eq.(2) is given by

$\displaystyle
J_{12}(\theta)=J_{12}^{0}\bigg{(}1-3\mathrm{cos}^{2}(\frac{\pi}{2}-\theta)\bigg{)}.$
(3)

Considering the schematic of reaction center as illustrated in Figure 1(a),
initially both donor chromophores are quite optically active, synergistically
facilitating the transition of excited electrons to the acceptor A. Dipole-
dipole interaction occurs because of the electrostatic attraction between the
positive and negative charges of the dipoles, impacting the electric field
experienced by them, which in turn affects energy levels, transition
probabilities, and other characteristics. The most important phenomena among
these photophysical properties is the orbital overlap between adjacent donor
molecules, resulting formation of new optically excitable states. Stable
delocalized excited states, known as exciton states, are commonly observed in
pigment-protein complexes [29]. These states play a crucial role in producing
quantum interference effects that boost the photocurrent and the power of our
QHE. In the presence of dipolar excitonic coupling $J_{12}$, the new donor
eigenstates are formed as

$\displaystyle|x_{1}\rangle=\frac{1}{\sqrt{2}}(|a_{1}\rangle+|a_{2}\rangle),$
(4)

which is symmetric/bright combinations of the uncoupled donor states and

$\displaystyle|x_{2}\rangle=\frac{1}{\sqrt{2}}(|a_{1}\rangle-|a_{2}\rangle),$
(5)

which is antisymmetric/dark combinations of the uncoupled donor states. Also
the corresponding energy eigenvalues are

$\displaystyle
E_{x_{1}}=\frac{E_{1}+E_{2}}{2}+\sqrt{\frac{(E_{1}-E_{2})^{2}}{4}+J_{12}^{2}(\theta)},$
(6)

and

$\displaystyle
E_{x_{2}}=\frac{E_{1}+E_{2}}{2}-\sqrt{\frac{(E_{1}-E_{2})^{2}}{4}+J_{12}^{2}(\theta)},$
(7)

where $E_{1}$ and $E_{2}$ represent the energies of the uncoupled donor states
$|a_{1}\rangle$ and $|a_{2}\rangle$, respectively. Due to the angle dependence
of the mentioned energy eigenvalues, it is obvious that different aggregation
structures will lead to distinct excited states and thus induce different
luminescent properties. Herein, we assume two coupled donors ($D_{1}$ and
$D_{2}$) are identical and degenerate, so the relationship between the angle
and the Davidov energy splitting between the symmetric and antisymmetric
states is given by $\bigtriangleup E(\theta)=2|J_{12}(\theta)|$. In the
following, we want to discuss about the aggregation structures named the
H-aggregate and the J-aggregate structure. In H-aggregate molecular
structures, the two transition dipole moments are aligned in head-to-head
manner with respect to vertical axis at $\theta=0$. Here, the symmetric state
is higher than the antisymmetric state (see Figure 1(b)), so the optical
transition is shifted to blue. On the other hand, in J-aggregate molecular
structures, the two transition dipole moments are aligned in head-to-tail
manner with respect to vertical axis at $\theta=\pi/2$ and the antisymmetric
state is higher than the symmetric state (see Figure 1(c)), so the optical
transition is changed to red. It also should be noted that the shift from the
H-aggregate to the J-aggregate occurred at the approximately critical angle
$\theta_{c}\thickapprox 35.26^{{}^{\circ}}$ with respect to vertical axis
[20]. It is worth mentioning that classically the total dipole moment is
always $2|\mu|$ when two dipole moments point to the same direction. But, the
dipole moment of $|x_{1}\rangle$ is strengthened through constructive
interference of the individual transition dipole matrix elements,
$\mu_{x_{1}}=(\mu_{1}+\mu_{2})/\sqrt{2}=\sqrt{2}|\mu|$, whereas the dipole
moment of $|x_{2}\rangle$ is annihilated due to destructive interference
$\mu_{x_{2}}=0$. Therefore, the symmetric combination describes an optically
active bright state with the optical transition rate
$\gamma_{h}\propto|\mu_{x_{1}}|^{2}=2|\mu|^{2}$, which is doubled in
comparison with an uncoupled donor case, while the antisymmetric combination
describes an optically forbidden dark state.

Unlike previous studies in which the starting point of the cyclic light-
emitting engine was considered in only ground state, in our new scheme, the
initial state of our model is assumed to be in the following desired state
given by

$\displaystyle|\psi\rangle=\alpha|x_{1}\rangle+\beta|x_{2}\rangle.$ (8)

Amazingly, we will find that under realistic constraints, different initial
state effects surpass the performance of previously investigated systems and
introduce an effective novel band gap engineering mechanism in order to
harness quantum effects and coherent superposition advantages. Next, the
excited electrons can be transferred to the acceptor molecule through
electronic coupling and the emission of phonons, following the process
described in [15]. Subsequently, the excited electrons are utilized to perform
work, resulting the charge-separated state $|\alpha\rangle$ decaying to the
sub-stable state $|\beta\rangle$. Also, the recombination between the acceptor
and the donor is also considered with a decay rate of $\chi\Gamma$, where
$\chi$ is a dimensionless fraction. Ultimately, the state $|\beta\rangle$
undergoes decay, returning via $\Gamma_{c}$ to the charge neutral ground state
$|b\rangle$ and complete the cycle.

We will now develop a kinetic scheme to describe the time evolution of the
average occupations associated with this specific level structure. As usual,
we use the standard Born-Markov approximations: a weak interaction between an
open quantum system and the environment and the extremely short correlation
time of the environment, that is, no memory effect. So, we can obtain the
master equation for the density matrix $\rho(t)$ of the photocell as follows

$\displaystyle\dfrac{\partial\rho}{\partial
t}=(\dfrac{i}{\hbar})[\rho,\mathcal{H}]+\mathcal{L}^{h}(\rho)+\mathcal{L}^{c}(\rho)+\mathcal{L}^{D}(\rho).$
(9)

The Hamiltonian of the photocell is given by
$\mathcal{H}=\mathcal{H}_{D}+\mathcal{H}_{A}$ where
$\mathcal{H}_{A}=\hbar\omega_{\alpha}\sigma_{\alpha}^{+}\sigma_{\alpha}^{-}+\hbar\omega_{\beta}\sigma_{\beta}^{+}\sigma_{\beta}^{-}$
is the Hamiltonian of the accepter and
$\sigma_{\alpha}^{+}=(\sigma_{\alpha}^{-})^{\dagger}=|\alpha\rangle\langle b|$
and $\sigma_{\beta}^{+}=(\sigma_{\beta}^{-})^{\dagger}=|\beta\rangle\langle
b|$ are the related raising (lowering) operators. The Lindblad operator,
denoted as $\mathcal{L}^{h}(\rho)$, encapsulates the interaction between the
system and its surrounding hot bath, it is given by

$\displaystyle\mathcal{L}^{h}(\rho)=(\dfrac{\gamma_{h}}{2})(n_{h}+1)\bigg{[}2\zeta^{-}\rho\zeta^{+}-\zeta^{+}\zeta^{-}\rho-\rho\zeta^{+}\zeta^{-}\bigg{]}+(\dfrac{\gamma_{h}}{2})n_{h}\bigg{[}2\zeta^{+}\rho\zeta^{-}-\zeta^{-}\zeta^{+}\rho-\rho\zeta^{-}\zeta^{+}\bigg{]},$
(10)

where $\zeta^{+}=(\zeta^{-})^{\dagger}=|x_{1}\rangle\langle b|$ is the
corresponding raising (lowering) operator. In addition, $\gamma_{h}$ is the
transition rate between ground state $|b\rangle$ and the bright state
$|x_{1}\rangle$ and $n_{h}$ is the average thermal occupations at $T_{h}$
according to the Planck distribution [30, 31, 32]. Likewise, the Lindblad
operator $\mathcal{L}^{c}(\rho)$, arising from the interaction with the cold
bath can be written as follows

$\displaystyle\mathcal{L}^{c}(\rho)=\mathcal{L}^{c_{1}}(\rho)+\mathcal{L}^{c_{2}}(\rho)+\mathcal{L}^{c_{3}}(\rho),$
(11)

where $\mathcal{L}^{c_{1}}(\rho)$ has the same form to Eq.(10) except
replacing $\gamma_{h}\rightarrow\gamma_{x}$, $n_{h}\rightarrow n_{x}$ and
$\zeta^{+}\rightarrow\xi^{+}$ with
$\xi^{+}=(\xi^{-})^{\dagger}=|x_{1}\rangle\langle x_{2}|$. Similarly, the
Lindblad operator $\mathcal{L}^{c_{2}}(\rho)$ ($\mathcal{L}^{c_{3}}(\rho)$)
has the same form to Eq.(10) except replacing
$\gamma_{h}\rightarrow\gamma_{c}$ ($\gamma_{h}\rightarrow\Gamma_{c}$),
$n_{h}\rightarrow n_{c}$ ($n_{h}\rightarrow N_{c}$) and
$\zeta^{+}\rightarrow\eta^{+}$ ($\zeta^{+}\rightarrow\tau^{+}$) where
$\eta^{+}=(\eta^{-})^{\dagger}=|x_{2}\rangle\langle\alpha|$
($\tau^{+}=(\tau^{-})^{\dagger}=|\beta\rangle\langle b|$) are the jump
operators. Correspondingly, $n_{x}$, $n_{c}$ and $N_{c}$ are the average
thermal occupations at $T_{c}$ with energies $\Delta E=E_{x_{1}}-E_{x_{2}}$,
$\Delta E=E_{x_{2}}-E_{\alpha}$ and $\Delta E=E_{\beta}-E_{b}$, respectively.
The average thermal occupations at energies $\Delta E$ are given as follows

$\displaystyle n=\dfrac{1}{e^{\Delta E/k_{B}T_{c}}-1},$ (12)

where $k_{B}$ stands for the Boltzmann constant. Ultimately, the Lindblad
operator $\mathcal{L}^{D}(\rho)$ is expressed in a more detailed form as

$\displaystyle\mathcal{L}^{D}(\rho)=(\dfrac{\Gamma}{2})\bigg{[}2\lambda_{1}^{-}\rho\lambda_{1}^{+}-\lambda_{1}^{+}\lambda_{1}^{-}\rho-\rho\lambda_{1}^{+}\lambda_{1}^{-}\bigg{]}+(\dfrac{\chi\Gamma}{2})\bigg{[}2\lambda_{2}^{+}\rho\lambda_{2}^{-}-\lambda_{2}^{+}\lambda_{2}^{-}\rho-\rho\lambda_{2}^{+}\lambda_{2}^{-}\bigg{]}.$
(13)

In the mentioned equation,
$\lambda_{1}^{+}=(\lambda_{1}^{-})^{\dagger}=|\alpha\rangle\langle\beta|$ and
$\lambda_{2}^{+}=(\lambda_{2}^{-})^{\dagger}=|\alpha\rangle\langle b|$ are
jump operators. The transition rate $\Gamma$ describes transitions from the
state $|\alpha\rangle$ to the state $|\beta\rangle$. The parameter $\chi$
characterizes the recombination of charge carriers (electrons and holes) at
the acceptor-donor interface effecting the current generation in a photocell.
By assuming $\chi=0$ we analyze the scenario where no recombination occurs.
This consideration aims to achieve the maximum power output from the
photocell. The kinetics of light-matter interactions, following the Pauli
master equation, guarantee completely positive populations. It is assumed that
our proposed scheme functions as a quantum heat engine while being in thermal
equilibrium with both hot and cold baths at the same time.

##  III. RESULTS AND DISCUSSION

Before the main discussion, we need to solve numerically the Pauli master
equation Eq.(9) and obtain the transient and steady states. In our
calculation, we use the parameters listed in Table 1.

Table 1. The used parameters.
---
$E_{a_{1}}-E_{b}$ | 1.8 (eV)
$E_{a_{2}}-E_{b}$ | 1.8 (eV)
$E_{a_{1}}-E_{\alpha}$ | 0.2 (eV)
$E_{a_{1}}-E_{\alpha}$ | 0.2 (eV)
$E_{\beta}-E_{b}$ | 0.2 (eV)
$\gamma_{h}=2\gamma_{1h}=2\gamma_{2h}$ | $1.24\times 10^{-6}$
$\gamma_{c}=2\gamma_{1c}=2\gamma_{2c}$ | $12\times 10^{-3}$
$\gamma_{x}$ | $25\times 10^{-3}$
$\Gamma_{c}$ | 0.0248
$J_{12}^{0}$ | 0.015

By incorporating the idea of photochemical current and voltage, we assign an
effective current and voltage to the reaction center. The resistance of an
external load is described by the electron decay rate $\Gamma$ from cathode to
anode. The current flows from $|\alpha\rangle$ to $|\beta\rangle$ is expressed
as $I=e\Gamma P_{\alpha}$, where $e$ represents the fundamental charge of an
electron at steady state and the probability to find the system in state
$i(i=x_{1},x_{2},\alpha,\beta,b)$ is denoted as $P_{i}=\rho_{i,i}$, where
$\rho_{i,i}$ is the diagonal elements of the density matrix. It should be
noted that the value of $\Gamma$ ranges from $\Gamma=0$ (the open circuit
regime) to large $\Gamma$ (the short circuit regime). The voltage across the
solar cell is determined by the chemical potential variance between the
mentioned two loads, $eV\equiv
e(V_{\beta}-V_{\alpha})$=$\mu_{\alpha}-\mu_{\beta}$. By employing the
Boltzmann distributions for levels $|\alpha\rangle$ and $|\beta\rangle$, where
$P_{\alpha}=e^{-(E_{\alpha}-\mu_{\alpha})/K_{B}T_{c}}$ and
$P_{\beta}=e^{-(E_{\beta}-\mu_{\beta})/K_{B}T_{c}}$ the solar cell’s voltage
is described in relation to energy levels and populations [17, 19]

$\displaystyle
eV=E_{\alpha}-E_{\beta}+K_{B}T_{c}\ln\dfrac{P_{\alpha}}{P_{\beta}}.$ (14)

When sunlight is absent, the system reaches thermal equilibrium with the
phonon bath temperature, resulting in a vanishing voltage. Consequently, $V$
serves as a metric for quantifying the deviation from the thermal state
characterized by temperature $T_{c}$. The rapid transfer rate of excited
electrons to the acceptor results in a current proportional to the product of
the electron charge and the generation rate, typically on the order of
microamperes. It should be noted that the current and voltage are determined
by analyzing the steady-state solutions of the master equation where we
calculate the populations $\rho_{\alpha,\alpha}$ and $\rho_{\beta,\beta}$ at
sufficiently long times. Now, the power of the quantum photocell is given by

$\displaystyle P_{\mathrm{out}}=VI.$ (15)

In this paper, in order to reach the desired power enhancement we integrate
aggregation optimization, system preparation in different initial delocalized
states, dark state protection and bandgap engineering. We discuss how
aggregation becomes important and initializing the quantum photocell in a
delocalized state gives rise to the big improvement of system performance
compared with the previous studies. In this regard, we define the enhanced
power output induced by the mentioned quantum effects as follows

$\displaystyle
PE=(\frac{P^{\mathrm{coupled}}_{\mathrm{out}}-P^{\mathrm{uocoupled}}_{\mathrm{out}}}{P^{\mathrm{uocoupled}}_{out}})\times
100,$ (16)

where $P^{\mathrm{coupled}}$ and $P^{\mathrm{uocoupled}}$ being the power
output in the coupled ($J_{12}\neq 0$) and uncoupled ($J_{12}=0$) case,
respectively.

Figure 3 shows the power generated by the quantum photocell as a function of
the coefficient $\beta$ in excitonically coupled $(J_{12}\neq 0)$ case for the
H-aggregate configuration. Extracting power has been evaluated at room
temperature $T_{c}=300$K. Figure 3(a) is plotted for the operation near the
short circuit regime $\Gamma=0.1$eV and Figure 3(b) is plotted for the open
circuit regime $\Gamma=0.001$eV. According to the Figure 3(a), when the
initial state is prepared in the ground state $|b\rangle$, the power output
always is higher compared to the general delocalized initial state
$|\psi\rangle$. Meanwhile, this result is completely reversed in the case of
$\Gamma=0.001$eV. The maximum power output
$P_{\mathrm{out}}^{\mathrm{max}}=4334.3\mu W$ of the studied photocell is
obtained for $\beta=0$ representing the initial state $|x_{1}\rangle$ in the
H-aggregation. Also, the production of output power values in descending order
according to the preparation of the specific initial states in H-aggregation
is: $P_{x_{1}}>P_{a_{1}}=P_{a_{2}}>P_{x_{2}}$. In contrast to the previous
claim, considering the H-aggregation type for $\Gamma=0.001eV$, when the
system is initialized in the general delocalized state $|\psi\rangle$ the
resulting power output is greater than that obtained from the ground state
$|b\rangle$. In contrast to the previous case, the power output of the
photocell reaches its minimum value for the initial state $|x_{1}\rangle$ and
the maximum power corresponds to the initial state $|x_{2}\rangle$.

Figure 4 presents the generated power as a function of the coefficient $\beta$
in excitonically coupled $(J_{12}\neq 0)$ case for the J-aggregate
configuration. As depicted in Figure 4(a,b), the analysis examines two
distinct values near the short circuit regime ($\Gamma=0.1$eV) and near the
open circuit regime ($\Gamma=0.001$eV), respectively. According to Figure
4(a), the power output is greater when the initial state is prepared in the
ground state $|b\rangle$ compared to the prepared initial state in the general
delocalized state $|\psi\rangle$. Conversely, this outcome is entirely
reversed in the case with $\Gamma=0.001$ eV. The minimum power output obtained
using the discussed quantum effects across various configurations and by
varying the $\Gamma$ parameter is $P_{\mathrm{out}}^{\mathrm{min}}=78.34\mu
W$, occurring with the initial state $|x_{2}\rangle$ in the J-aggregation and
for $\Gamma=0.001$eV. As observed in Figure 4(a), the descending order of
output power values based on the specific initial states preparation in the
J-aggregate configuration is as follows:
$P_{x_{2}}>P_{a_{1}}=P_{a_{2}}>P_{x_{1}}$. Similar to Figure 3(b) when the
system is initialized in the general delocalized state $|\psi\rangle$ the
resulting power output is always greater than that obtained from the ground
state $|b\rangle$. In this scenario, the power output of the photocell reaches
its minimum for the initial state $|x_{2}\rangle$. Conversely, the maximum
power output corresponds to the initial state $|x_{1}\rangle$ at $\beta=0$,
demonstrating the peak performance.

The angle dependence of the power enhancement $\mathrm{PE}$ near the short and
open circuit regimes is illustrated in Figure 5. As shown in Figure 5(a), it
is obvious that the power enhancement is largely independent of the system’s
initial state, as the graphs are nearly identical and overlaped. In fact, the
power enhancement does not significantly depend on the initial state and
starting point of the system. Moreover, the enhancement of power is affected
by the aggregation of dipole moments, with the optimal condition occurring at
the H-aggregate or when $\theta=0$. As $\theta$ increases, the power
enhancement decreases, hitting its minimum value in the J-aggregate state. At
the critical angle $\theta_{c}$ (transition from the H-aggregate phase to the
J-aggregate), $J_{12}$ becomes zero, indicating the possibility of an
uncoupled regime. Therefore, the uncoupled power is equal to the power of the
coupled regime, resulting in zero power enhancement. In the following, we
investigate the angle dependency of the $\mathrm{PE}$ with the used parameter
$\Gamma=0.001$eV, as depicted in Figure 5(b). In this scenario, the
$\mathrm{PE}$ interestingly exhibits a significant reliance on the initial
state of the system. For the H-aggregate case, the $\mathrm{PE}$ demonstrates
positive values and it follows a ascending trend for the different initial
$\beta=(0,2^{-1/2},1)$ with approximately values of $17.64\%$, $26.75\%$ and
$35.87\%$, respectively. It drops zero at magic angle $\theta=\theta_{c}$ and
after crossing the critical angle, the $\mathrm{PE}$ takes negative values and
ultimately it approximately reaches $-42.19\%$, $-64.64\%$ and $-87.08\%$ in
the J aggregate configuration of dipoles. Overall, we conclude that when
$\theta<\theta_{c}$, the power of coupled case in higher than the power of
uncoupled case and when $\theta>\theta_{c}$, the power of coupled case in
lower than the power of uncoupled case.

To give a more clear physical reason for our results in Figure 5, we emphasize
that for the H-aggregate, the optimal state occurs at $\theta=0$, resulting in
maximum power enhancement due to constructive dipole coupling. As $\theta$
increases from this optimal alignment, the power enhancement decreases,
reflecting a transition from constructive to destructive interference of
dipole interactions. Conversely, in the J-aggregates, the power enhancement
reaches its lowest value. This is due to the increased destructive
interference as the dipoles misaligned. Moreover, it is worth noting that for
the H-aggregate case ($\theta<\theta_{c}$), the bright state $|x_{1}\rangle$
lies at a higher energy than the dark state $|x_{2}\rangle$ and for the
J-aggregate case ($\theta>\theta_{c}$), the bright state $|x_{1}\rangle$ lies
at a lower energy than the dark state $|x_{2}\rangle$. Also, by comparing the
energy splitting between the bright and dark states for the H and
J-aggregates, we have $\Delta E^{\mathrm{J-aggregate}}=2\Delta
E^{\mathrm{H-aggregate}}$. Consequently, the transition from the donor to the
acceptor for the J-aggregate case is very low and the current enhancement is
negative because an electron in the bright state jumps to the dark state
through only the stimulated absorption $\gamma_{x}n_{x}$ of thermal phonons
for the J-aggregate case while this transition occurs through the stimulated
and spontaneous emissions of thermal phonons $\gamma_{x}(1+n_{x})$ for the
H-aggregate case.

##  IV. CONCLUSIONS

In summary, the importance of the initial-state dependence for a photovoltaic
cell efficiency under different aggregation and realistic constraints has been
investigated near the short and open circuit regimes. We have considered
dipole molecules of solar materials producing delocalized states and amazingly
found that the power output with the initially prepared general delocalized
excited state can surpass the performance of previously studied systems with
the initially prepared ground state near the open circuit regime and in the
both H and J aggregation. In addition, our findings indicated that our
considered photocell can result in a notable enhancement of power output by
about $35.87\%$ in a coherent coupled dipole system compared with the
uncoupled dipoles in the H aggregate and near the open circuit regime. This
study provides valuable insights into the improvement in PV cells conversion
efficiency offering a foundation for future investigations and practical
applications in light-harvesting industry.

## Data availability

The data of the present study are available from the corresponding author upon
a reasonable request.

## Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * [1] Tokihiro, Tetsuji, Yoshihiro Manabe, and Eiichi Hanamura. ”Superradiance of Frenkel excitons in linear systems.” Physical Review B 47.4 (1993): 2019.
  * [2] Calhoun, Tessa R., et al. ”Quantum coherence enabled determination of the energy landscape in light-harvesting complex II.” The Journal of Physical Chemistry B 113.51 (2009): 16291-16295.
  * [3] Abramavicius, Darius, Benoit Palmieri, and Shaul Mukamel. ”Extracting single and two-exciton couplings in photosynthetic complexes by coherent two-dimensional electronic spectra.” Chemical physics 357.1-3 (2009): 79-84.
  * [4] Panitchayangkoon, Gitt, et al. ”Long-lived quantum coherence in photosynthetic complexes at physiological temperature.” Proceedings of the National Academy of Sciences 107.29 (2010): 12766-12770.
  * [5] Harel, Elad, Andrew F. Fidler, and Gregory S. Engel. ”Real-time mapping of electronic structure with single-shot two-dimensional electronic spectroscopy.” Proceedings of the National Academy of Sciences 107.38 (2010): 16444-16447.
  * [6] Hayes, Dugan, Graham B. Griffin, and Gregory S. Engel. ”Engineering coherence among excited states in synthetic heterodimer systems.” Science 340.6139 (2013): 1431-1434.
  * [7] Romero, Elisabet, et al. ”Quantum coherence in photosynthesis for efficient solar-energy conversion.” Nature physics 10.9 (2014): 676-682.
  * [8] Mohseni, Masoud, et al. ”Environment-assisted quantum walks in energy transfer of photosynthetic complexes.” J. Chem. Phys 129.17 (2008): 174106.
  * [9] Plenio, Martin B., and Susana F. Huelga. ”Dephasing-assisted transport: quantum networks and biomolecules.” New Journal of Physics 10.11 (2008): 113019.
  * [10] Rebentrost, Patrick, et al. ”Environment-assisted quantum transport.” New Journal of Physics 11.3 (2009): 033003.
  * [11] Zhu, Jing, et al. ”Modified scaled hierarchical equation of motion approach for the study of quantum coherence in photosynthetic complexes.” The Journal of Physical Chemistry B 115.6 (2011): 1531-1537.
  * [12] Yeh, Shu-Hao, Jing Zhu, and Sabre Kais. ”Population and coherence dynamics in light harvesting complex II (LH2).” The Journal of Chemical Physics 137.8 (2012).
  * [13] Fenna, R. E., and B. W. Matthews. ”Chlorophyll arrangement in a bacteriochlorophyll protein from Chlorobium limicola.” Nature 258.5536 (1975): 573-577.
  * [14] Blankenship, Robert E., et al. ”Comparing photosynthetic and photovoltaic efficiencies and recognizing the potential for improvement.” science 332.6031 (2011): 805-809.
  * [15] Dorfman, Konstantin E., et al. ”Photosynthetic reaction center as a quantum heat engine.” Proceedings of the National Academy of Sciences 110.8 (2013): 2746-2751.
  * [16] Scully, Marlan O., et al. ”Quantum heat engine power can be increased by noise-induced coherence.” Proceedings of the National Academy of Sciences 108.37 (2011): 15097-15100.
  * [17] Scully, Marlan O. ”Quantum photocell: Using quantum coherence to reduce radiative recombination and increase efficiency.” Physical review letters 104.20 (2010): 207701.
  * [18] Svidzinsky, Anatoly A., Konstantin E. Dorfman, and Marlan O. Scully. ”Enhancing photovoltaic power by Fano-induced coherence.” Physical Review A 84.5 (2011): 053818.
  * [19] Creatore, Celestino, et al. ”Efficient biologically inspired photocell enhanced by delocalized quantum states.” Physical review letters 111.25 (2013): 253601.
  * [20] Zhang, Yiteng, et al. ”Delocalized quantum states enhance photocell efficiency.” Physical chemistry chemical physics 17.8 (2015): 5743-5750.
  * [21] Fruchtman, Amir, et al. ”Photocell optimization using dark state protection.” Physical review letters 117.20 (2016): 203603.
  * [22] Friend, Richard H., et al. ”Electroluminescence in conjugated polymers.” Nature 397.6715 (1999): 121-128.
  * [23] Uoyama, Hiroki, et al. ”Highly efficient organic light-emitting diodes from delayed fluorescence.” Nature 492.7428 (2012): 234-238.
  * [24] Samuel, Ifor David Williams, and Graham Alexander Turnbull. ”Organic semiconductor lasers.” Chemical reviews 107.4 (2007): 1272-1295.
  * [25] Tessler, Nir. ”Lasers based on semiconducting organic materials.” Advanced Materials 11.5 (1999): 363-370.
  * [26] Hu, Wenping, et al. Organic optoelectronics. John Wiley & Sons, 2012.
  * [27] Dong, Huanli, et al. ”Organic photoresponse materials and devices.” Chemical Society Reviews 41.5 (2012): 1754-1808.
  * [28] Sirringhaus, Henning. ”25th anniversary article: organic field‐effect transistors: the path beyond amorphous silicon.” Advanced materials 26.9 (2014): 1319-1335.
  * [29] Van Amerongen, Herbert, and Rienk Van Grondelle. Photosynthetic excitons. World Scientific, 2000.
  * [30] Carmichael, H. J. ”An open systems approach to quantum optics Springer.” Lecture notes in physics 18 (1993).
  * [31] Breuer, Heinz-Peter, and Francesco Petruccione. The theory of open quantum systems. OUP Oxford, 2002.
  * [32] Mandel, Leonard, and Emil Wolf. Optical coherence and quantum optics. Cambridge university press, 1995.

$\mathbf{(a)}$ $\mathbf{(b)}$

$\mathbf{(c)}$ $\mathbf{(d)}$

$\mathbf{(e)}$ $\mathbf{(f)}$

Figure 1: Schematic of our PV cell with the energy level diagrams. In (a),
uncoupled case, both donors $D_{1}$ and $D_{2}$ exhibit optical activity and
synergistically facilitate the transfer of excited electrons to the acceptor
A. The red and blue shaded/pale regions encircling the molecules illustrate
the molecular orbitals, which depict the spatial electron density
distribution. In (c) and (e), coupled cases, the interaction between $D_{1}$
and $D_{2}$ leads to the formation of a coupled system which generates new
eigenstates superposed of the uncoupled donor states $|a_{1}\rangle$ and
$|a_{2}\rangle$. (b), (d) and (f) show the level structure and electron paths
through bright and dark states of uncoupled, H-aggregate and J-aggregate,
respectively.

Figure 2: The alignment of two parallel dipole moments $\mu_{1}$ and $\mu_{2}$
is depicted as yellow arrows. Here, the angle $\theta$ is measured with
respect to the vertical axis.

(a)

(b)

Figure 3: Power generated in PV cell as a function of the coefficient $\beta$
at room temprature when system reaches steady state operation in the presence
of coherent dipole coupling ($J_{12}\neq 0$) with the H-aggregate
configuration for (a) $\Gamma=0.1$eV and (b) $\Gamma=0.001$eV. The red dashed
line is a benchmark line denoted for comparison of the generated power in all
delocalized states with the initially prepared ground state $|b\rangle$.

(a)

(b)

Figure 4: Power generated in PV cell as a function of the coefficient $\beta$
at room temprature when system reaches steady state operation in the presence
of coherent dipole coupling ($J_{12}\neq 0$) with the J-aggregate
configuration for (a) $\Gamma=0.1$eV and (b) $\Gamma=0.001$eV. The red dashed
line is a benchmark line denoted for comparison of the generated power in all
delocalized states with the initially prepared ground state $|b\rangle$.

(a) (b)

Figure 5: Power enhancement as a function of the angle $\theta$ with respect
to the vertical axis for (a) $\Gamma=0.1$eV and (b) $\Gamma=0.001$eV. Blue
circle line denotes power proportional to $\beta=0$, green square line denotes
$\beta=2^{-1/2}$ and red line denotes $\beta=1$ representing the set up
prepared in initial states $|x_{1}\rangle$, $|a_{1}\rangle$ and
$|x_{2}\rangle$, respectively. The black arrow indicates the specific angle
known as the magic angle $\theta_{c}\approx 35.26^{\circ}$.