# Winning the CityLearn Challenge: Adaptive Optimization
with Evolutionary Search under Trajectory-based Guidance

Vanshaj Khattar and Ming Jin

###### Abstract

> Modern power systems will have to face difficult challenges in the years to
> come: frequent blackouts in urban areas caused by high power demand peaks,
> grid instability exacerbated by intermittent renewable generation, and
> global climate change amplified by rising carbon emissions. While current
> practices are growingly inadequate, the path to widespread adoption of
> artificial intelligence (AI) methods is hindered by missing aspects of
> trustworthiness. The CityLearn Challenge is an exemplary opportunity for
> researchers from multiple disciplines to investigate the potential of AI to
> tackle these pressing issues in the energy domain, collectively modeled as a
> reinforcement learning (RL) task. Multiple real-world challenges faced by
> contemporary RL techniques are embodied in the problem formulation. In this
> paper, we present a novel method using the solution function of optimization
> as policies to compute actions for sequential decision-making, while notably
> adapting the parameters of the optimization model from online observations.
> Algorithmically, this is achieved by an evolutionary algorithm under a novel
> trajectory-based guidance scheme. Formally, the global convergence property
> is established. Our agent ranked first in the latest 2021 CityLearn
> Challenge, being able to achieve superior performance in almost all metrics
> while maintaining some key aspects of interpretability.

##  1 Introduction

Rapid urbanization in the past decades has led to a substantial increase in
energy use that puts stress on the grid assets, while the integration of
additional renewable generation and energy storage at the distribution level
introduces both opportunities and new challenges (Rolnick et al. 2022). The
cornerstone of addressing emerging issues is the deployment of proper control
and coordination strategies, which have a potential impact on enhancing energy
flexibility and resilience in the face of a surge in climate-induced demand
(as already seen in places like California, where rolling blackouts are
increasingly frequent during the summer) (Vazquez-Canteli et al. 2020).
Current industry practice is heavily based on optimization models, such as
energy dispatch and unit commitment, where parameters (e.g., technological and
physical constraints) are fixed throughout the lifecycle; however, such an
approach is increasingly confronted by environmental uncertainty, renewable
generation stochasticity, and the ever-increasing complexity of the
distribution grid (Abedi, Gaudard, and Romerio 2019). On the other hand, there
has been a surge in machine learning research, notably RL, because it allows
the agent to act without the need to access the true model—a feature of
particular interest for large-scale, complex systems, where it is not cost-
effective to develop models of such high fidelity. Despite recent progress,
real-world RL is still in its infancy (Dulac-Arnold et al. 2021).

Against this backdrop, the CityLearn Challenge aims to spur RL solutions for
the control of modern energy systems by providing a set of benchmarks for
urban energy management, load shaping, and demand response in a range of
climate zones (Vazquez-Canteli et al. 2020). The agent is tasked with
exploring and exploiting the best control strategy for energy storage
distributed in a community of buildings. Performance is evaluated against
standard metrics such as ramping costs, peak demands, and carbon emissions.
The CityLearn encapsulates 4 of the 9 real-world RL challenges identified by
(Dulac-Arnold et al. 2021), including _1)_ the ability to learn on live
systems from limited samples—there is no training period; _2)_ dealing with
system constraints that should never or rarely be violated—there are strict
balancing equations for electricity, heating, and cooling energy; _3)_ the
ability to provide quick action—there is a strict time limit for completing
the 4-year evaluation on Google’s Colab; and _4)_ providing system operators
with explainable policies—a necessity to facilitate real-world adoption and
deployment.

Figure 1: Overview. (A) Agent trajectory of observed states (square) and
actual (green triangle)/planned (grey triangle) actions. At each step, the
agent solves optimization (2) to plan ahead but only executes the most
immediate action; $\bar{a}_{t:t^{\prime}}$ represents the action planned at
time $t$ for time $t^{\prime}\geq t$. (B) Illustration of the system,
including interaction with the environment and adaptation of optimization
parameters. Each iteration of $k$ consists of evaluating $N_{k}$ (three in
this case) policies sampled from the candidate distribution. The observed
rewards $\\{R_{j}^{k}\\}_{j\in[N_{k}]}$ and trajectories
$\\{s_{t}^{j,k},a_{t}^{j,k}\\}_{j\in[N_{k}]}$ are stored in some buffers,
which are then used to compute weights $\\{\nu_{j}^{k}\\}_{j\in[N_{k}]}$ and
guidance signals $\\{\varrho_{j}^{k}\\}_{j\in[N_{k}]}$ to update the
distribution as in (5). (C) The two-level structure indicates the
correspondence between an optimization parameter $\zeta$ and the solution
function $\pi_{\zeta}$ (used as policy). The upper-level objective is the
expected episodic rewards (1) and the lower-level objective is the objective
function of (2), the extremum of which is the policy action (see (3)). The
guidance signal $\varrho$ is computed based on the trajectory of each
candidate and applied on top of the original parameter during the update.

In this paper, we describe our winning solution for the 2021 CityLearn
Challenge based on the idea broadly categorized as _adaptive optimization_.
Indeed, optimization (especially convex optimization) has become the de facto
standard in industrial systems with profound theoretical foundations, and
various formulations for control and planning applications (Boyd, Boyd, and
Vandenberghe 2004). Such approaches can easily encode domain-specific
constraints (in the form of nonlinear functions and variational inequalities)
and can gracefully handle problems with millions of decision variables
(Facchinei and Pang 2007). Although well established, optimization models,
once built, typically do not adapt to real-world conditions, rendering current
approaches rather “rigid.”

To address this fundamental limitation, we exploit that the solution of
optimization lies on a manifold implicitly defined by a general equation
(Dontchev and Rockafellar 2009). The crux of our idea is to shape this
manifold by adapting the parameters of the optimization model while extracting
insights from trajectory data to design guidance signals (see Fig. 1 for a
detailed illustration). Key differences between our method and well-
established optimization techniques (e.g., stochastic optimization (Powell
2020), bi-level optimization (Dempe and Zemkoho 2020)) include:

  1. _1)_

we only allow access to the environment through interactive samples (reward,
states, etc.) but not the true dynamics or reward function;

  2. _2)_

we make full use of the control trajectory of a Markov decision process (MDP)
to obtain a guidance signal.

While _1)_ is similar to the RL setup, _2)_ can be viewed as an _augmentation
of zeroth-order search methods with insights extracted from control data_ ,
which, to the best of our knowledge, is the _first of its kind_. In general,
zeroth-order algorithms, such as simultaneous perturbation (Spall 2005; Mania,
Guy, and Recht 2018), evolutionary algorithms (Salimans et al. 2017; Zhou, Yu,
and Qian 2019), and Bayesian optimization (Snoek, Larochelle, and Adams 2012;
Frazier 2018) are natural candidates for RL and easy to implement, but can
potentially suffer scalability problems (Ghadimi and Lan 2013). Nevertheless,
the parameters of an optimization model (i.e., variables to be learned)
usually have clear interpretations. Thus, we design a mechanism to leverage
domain knowledge for the design of appropriate guidance during the search for
these parameters; such guidance can be specified as a general function of
trajectory data (including observed states and actions of an MDP). The method
works well in an online environment without an extensive training period,
which is particularly advantageous in a real-world setting where an offline
environment for model training is usually not available. According to
independent evaluations, the proposed method achieved the highest performance
in the recent 2021 CityLearn Challenge. To demonstrate effectiveness against
existing techniques, we further validate the method by comparing it with a
range of baselines. Key contributions are as follows:

  * •

A framework of adaptive optimization for online control, winning the 1st place
in the 2021 CityLearn Challenge

  * •

A novel evolutionary search (ES) algorithm with a guidance function based on
state-action-trajectory data

  * •

Theoretical analysis of asymptotic convergence to global optima with noisy
function evaluations

  * •

Empirical comparison against a range of baselines

###  1.1 Related work

Optimal control and stochastic optimal control are well-known approaches to
sequential decision-making problems (Bertsekas 2019). Convex optimization is
another avenue (Agrawal et al. 2020). Most existing works assume a known
dynamic model of the system. Various large-scale stochastic programs have been
proposed in the literature to deal with future uncertainty (Prékopa 2013;
Powell 2020). The major drawback is the computational burden of rapidly
expanding scenario trees in multi-stage stochastic programming. Our method
relieves computation by using plug-in estimators, a.k.a., deterministic
approximation of future uncertainty within convex optimization.

Recently, RL has gained popularity in various domains (Chen et al. 2022;
Haydari and Yilmaz 2020; Nian, Liu, and Huang 2020). To contextualize the
present approach, we make a few remarks about the relation to model-based RL
(MBRL). In particular _implicit_ MBRL, where the entire procedure (e.g., model
learning and planning) is optimized for optimal policy computation (Moerland,
Broekens, and Jonker 2020). However, unlike existing works (e.g., (Tamar et
al. 2016; Karkus, Hsu, and Lee 2017; Racanière et al. 2017; Guez et al. 2018;
Schrittwieser et al. 2020) that build a model based on
(recurrent/convolutional) neural networks (NNs) with primary restrictions to
discrete state and action space, our method learns how to plan by solving
optimization and adapting its parameters; hence, it is amenable to a wide
range of applications with continuous states and actions. The present work is
closely related to (Ghadimi, Perkins, and Powell 2020; Agrawal et al. 2020),
which also use convex optimization as a policy class to handle uncertainty. In
particular, convex optimization control policies are learned in (Agrawal et
al. 2020) with implicit differentiation (Agrawal et al. 2019). We extend their
method to the RL setting and propose a novel ES algorithm for learning.

Differentiated from existing ES-based strategies (Szita and Lörincz 2006;
Salimans et al. 2017; Khadka and Tumer 2018; Gangwani and Peng 2018; Conti et
al. 2018) or zeroth-order optimization (Snoek, Larochelle, and Adams 2012; Liu
et al. 2020), our ES is augmented with a guidance function that depends on the
data collected as the policy interacts with the environment; hence, it can be
viewed as a type of _MDP-augmented ES_. The guidance mechanism is also
flexible enough to allow the effective incorporation of domain knowledge, as
demonstrated in CityLearn. We further provide theoretical justification for
this rather complex scheme.

##  2 Preliminaries

###  2.1 Problem setup

Consider an MDP $(\mathcal{S},\mathcal{A},\mathcal{P},r)$, where $\mathcal{S}$
is the (possibly infinite) state space, $\mathcal{A}$ is the set of actions,
$\mathcal{P}:\mathcal{S}\times\mathcal{A}\rightarrow\mathcal{M}(\mathcal{S})$
is the transition probability kernel with $\mathcal{M}(\mathcal{S})$ denoting
the set of all probability measures over $\mathcal{S}$, and $r(s,a)$ gives the
corresponding immediate reward (can be time-dependent). The goal in episodic
RL is to learn a policy $\pi:\mathcal{S}\rightarrow\mathcal{A}$ that maximizes
cumulative rewards over a finite time horizon:

$\displaystyle\underset{\pi\in\Pi}{\max}$
$\displaystyle\mathbb{E}\left[R(\pi)\right]{,}$ (1)

where $R(\pi)\coloneqq{\sum}_{t=0}^{T}\hskip
2.84544ptr_{t}\big{(}s_{t},\pi(s_{t})\big{)}$ is the episodic reward, with
$s_{t}\in\mathcal{S}\subseteq\mathbb{R}^{n_{s}}$ denoting the state at time
$t$ and $T$ as the horizon. The expectation is taken over initial state
distributions and transition dynamics (under deterministic policy $\pi$). We
require access to a random sample $R(\pi)$ of the episodic reward as well as
trajectory data $\\{(s_{t},a_{t})\\}_{t\in[T]}$ to be used to later compute
the guidance signal. Here, we use the shorthand $[T]=\\{1,...,T\\}$.

###  2.2 Canonical approaches and solution functions

The proposed method is based on canonical stochastic programming methods
(Powell 2020). For example, in multi-stage stochastic programming (Pflug and
Pichler 2014), the action at state $s_{t}$ is computed as

$\underset{a_{t}\in\mathcal{A}}{\operatorname*{arg\,max}}\Big{(}\tilde{r}_{t}(s_{t},a_{t})\\!+\\!\\!\\!\\!\underset{\\{a_{t^{\prime}}\\}_{t^{\prime}=t+1}^{T}}{\max}\\!\\!\\!\tilde{\mathbb{E}}\Big{[}\\!\\!{\sum}_{t^{\prime}=t+1}^{T}\\!\tilde{r}_{t^{\prime}}\big{(}s_{t^{\prime}},a_{t^{\prime}})\Big{|}s_{t},a_{t}\Big{]}\Big{)},$

where $\tilde{r}_{t}:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$ is the
surrogate reward function and $\tilde{\mathbb{E}}[\cdot]$ is the surrogate
expectation operator (e.g., model-based scenario trees) designed to
approximate the true environment. The above formulation can also be seen as
finding the solution to a Bellman equation in dynamic programming. The main
limitations, nevertheless, are the high computational cost of evaluating the
expectation operator and the potential model mismatch due to approximations.

A simpler yet more practically appealing method, widely adopted in today’s
industries, is to use deterministic approximations of the future and capture
the dependence of future states on prior decisions through constraints as part
of the lookahead model (Powell 2020): $\pi_{\zeta}(s_{t})$ is given by

$\displaystyle\underset{\bar{a}_{t}\in\mathcal{A}}{\operatorname*{arg\,max}}\max_{\\{\bar{s}_{t^{\prime}},\bar{a}_{t^{\prime}}\\}_{t^{\prime}=t+1}^{T}}\bar{R}(s_{t},\bar{a}_{t},\\{\bar{s}_{t^{\prime}},\bar{a}_{t^{\prime}}\\}_{t^{\prime}=t+1}^{T};\zeta)$
(2) $\displaystyle\hskip 0.56917pt\text{s. t.}\quad
g_{i}(s_{t},\bar{a}_{t},\\{\bar{s}_{t^{\prime}},\bar{a}_{t^{\prime}}\\}_{t^{\prime}=t+1}^{T};\zeta)\leq
0\hskip 4.83691pt;\;\;\hskip 2.56064pt\forall i\in\mathcal{I}$
$\displaystyle\hskip 1.13791pt\quad\;\quad
h_{i}(s_{t},\bar{a}_{t},\\{\bar{s}_{t^{\prime}},\bar{a}_{t^{\prime}}\\}_{t^{\prime}=t+1}^{T};\zeta)=0\hskip
5.40608pt;\;\;\;\forall i\in\mathcal{E},$

with the objective
$\bar{R}(s_{t},\bar{a}_{t},\\{\bar{s}_{t^{\prime}},\bar{a}_{t^{\prime}}\\}_{t^{\prime}=t+1}^{T};\zeta)$
defined as

$\bar{r}_{t}\big{(}{s_{t}},\bar{a}_{t};\zeta\big{)}+{\sum}_{t^{\prime}=t+1}^{T}\bar{r}_{t^{\prime}}\big{(}\bar{s}_{t^{\prime}},\bar{a}_{t^{\prime}};\zeta\big{)},$

where $\\{\bar{s}_{t^{\prime}}\\}_{t^{\prime}=t+1}^{T}$ and
$\\{\bar{a}_{t^{\prime}}\\}_{t^{\prime}=t}^{T}$ are optimization variables
corresponding to the planned states and actions, $\bar{r}_{t}$ is surrogate
reward, and the feasible set is defined by $\\{g_{i}\\}_{i\in\mathcal{I}}$ and
$\\{h_{i}\\}_{i\in\mathcal{E}}$. We denote the parameters of the objective
function and the constraints collectively by
$\zeta\in\mathcal{Z}\subset\mathbb{R}^{d}$. The dependencies of future states
on current and _planned_ states/actions are encoded as constraints in (2) as
part of the lookahead model. Many examples can be found in, e.g., (Borrelli,
Bemporad, and Morari 2017). We remark that some optimization parameters may be
provided by predictors based on the current state $s_{t}$, the parameter of
which is also collected by $\zeta$ in this case.

The policy $\pi_{\zeta}(s_{t})$, a.k.a., the solution function (Dontchev and
Rockafellar 2009), provides the action at the current state $s_{t}$ as the
optimal solution to (2). Since this function is generally set-valued (Dontchev
and Rockafellar 2009), we make the following assumption.

######  Assumption 1.

For each $\zeta\in\mathcal{Z}$ and $s_{t}\in\mathcal{S}$: _a)_ the objective
function in (2) is continuous, strictly convex, $g_{i}$ is continuous and
convex for each $i\in\mathcal{I}$, and $h_{i}$ is affine for each
$i\in\mathcal{E}$; _b)_ the feasible set of (2) is closed, absolutely bounded,
and has a nonempty interior.

The above assumption can be satisfied by imposing proper conditions on the
design of the surrogate model, i.e., objective and constraints in (2). Note
that in our approach, we make no convexity assumption about the true dynamics
or rewards of the environment, which can be seen as a blackbox. The convexity
condition is only stipulated for the surrogate model for computational
efficiency. Our goal is simply to learn the parameters of the optimization
model in order to have a good decision-making capability.111Perhaps
surprisingly, despite the fact that (2) is convex, the policy (given by the
solution function) can be highly nonconvex with high representational
capacity. In particular, the contemporary work by (Jin et al. 2023) shows a
“universal approximation” property of the solution functions of linear
programs (LPs).  An immediate consequence of the above assumption is that the
solution to (2) is unique; furthermore, it implies continuity with respect to
parameters.

######  Lemma 1.

The solution function $\pi_{\zeta}(s_{t})$ defined in (2) is continuous with
respect to parameter $\zeta$ for each $s_{t}\in\mathcal{S}$.

The proof is a direct application of the Berge maximum theorem (Berge 1997).
To conclude this section, let us make some comments on the construction of the
surrogate model in (2). By analogy to reward design (Prakash et al. 2020), the
objective function should be chosen to promote desirable behaviors. The set of
constraints introduces inductive bias on the transition dynamics of the
environment. It is beneficial, though oftentimes unlikely and non-essential,
that the surrogate model matches the functional forms of true reward or
dynamics, an idea shared in model-based RL (Moerland, Broekens, and Jonker
2020). It is, nevertheless, desirable to ensure the computational efficiency
of (2) to provide actions quickly—hence the choice of convex programs.

##  3 Policy adaptation with ES under guidance

The potential mismatch between the surrogate model and the real environment
and errors due to predictions may adversely affect the decision quality of
(2). Thereby, we aim to adapt the parameters of the surrogate model to shape
the solution function. The task of finding the optimal parameter within the
set of solution functions $\Pi=\\{\pi_{\zeta}:\zeta\in\mathcal{Z}\\}$ can be
compactly written as:

$\Upsilon\coloneqq\underset{\zeta\in\mathcal{Z}}{\operatorname*{arg\,max}}\hskip
1.42271pt\mathbb{E}\Bigg{[}{\sum}_{t=0}^{T}r_{t}\big{(}s_{t},\pi_{\zeta}(s_{t})\big{)}\Bigg{]},$
(3)

where $\Upsilon$ is the set of global optima for policy parameters. Note that
$\zeta$ is not part of the true reward (which remains unknown to the agent)
but only the parameters of the surrogate model that implicitly defines the
policy in (2). Since $\pi_{\zeta}(s_{t})$ is given by an optimization (2), (3)
can also be viewed as a bi-level problem (c.f., (Dempe and Zemkoho 2020)): the
outer level aims to learn parameters to maximize rewards, while the inner
level defines policy action as a solution to (2). The key challenge in solving
(3) as a bi-level problem, however, is that the outer level objective is not
analytically revealed (thus prohibiting any direct differentiation approach)
and can be nonconvex with respect to $\zeta$.

Input: Hyperparameters $\\{N_{k}\\}$, uniform distribution $\mu$, initial
point $z_{0}$

1: Initialize $P_{1}(d\zeta)\sim\exp(\|\zeta-z_{0}\|)\mu(d\zeta)$

2: for $k=1,2,\ldots$ do

3: Sample $N_{k}$ candidates from the distribution $p_{k}$:
$\zeta_{1}^{k},\zeta_{2}^{k},\cdots,\zeta_{N_{k}}^{k}\overset{\text{iid}}{\sim}p_{k}$

4: for $j=1,\ldots,N_{k}$ do

5: Deploy policy $\pi_{\zeta_{j}^{k}}$ for one episode and observe an episodic
reward $R_{j}^{k}\leftarrow R(\pi_{\zeta_{j}^{k}})$

6: Compute the guidance signal $\varrho_{j}^{k}$ by (4)

7: end for

8: Update the distribution $p_{k+1}$ for the next iteration according to (5).

9: end for

Algorithm 1 Evolutionary search under guidance

###  3.1 Guided evolutionary search

This section discusses the proposed ES algorithm (Algorithm 1), inspired by
the method of generations (Zhigljavsky 2012). At each iteration $k$, the
algorithm randomly samples a set of $N_{k}$ parameter candidates,
$\zeta_{1}^{k},\cdots,\zeta_{N_{k}}^{k}\overset{iid}{\sim}p_{k}$. For each
candidate $j\in[N_{k}]$, we evaluate the corresponding policy in the
environment and observe an episodic reward
$R_{j}^{k}\sim\mathcal{R}(\pi_{\zeta_{j}^{k}})$, where
$\mathcal{R}(\pi_{\zeta_{j}^{k}})$ denotes the distribution of episodic reward
for policy $\pi_{\zeta_{j}^{k}}$, as well as all the state and action pairs
$\\{s^{j,k}_{t},a^{j,k}_{t}\\}_{t\in[T]}$ in the past episode. Based on
trajectory data, we compute a guidance signal

$\varrho_{j}^{k}=\phi(\\{s^{j,k}_{t},a^{j,k}_{t}\\}_{t\in[T]}),$ (4)

where $\phi:(\mathcal{S}\times\mathcal{A})^{T}\rightarrow\mathcal{Z}^{\prime}$
is a function that may have complicated dependence on past states and actions
with $\mathcal{Z}^{\prime}$ as the range. The design of such a guidance
function is often based on domain knowledge (to be discussed later in
CityLearn). Then, we update the distribution for the next iteration as

$p_{k+1}(d\zeta)=\sum_{j=1}^{N_{k}}\nu_{j}^{k}Q_{k}(\zeta_{j}^{k},\varrho_{j}^{k},d\zeta),$
(5)

where

$\nu_{j}^{k}=\frac{\exp(R_{j}^{k})}{\sum_{j=1}^{N_{k}}\exp(R_{j}^{k})}$ (6)

are the weights obtained by taking the softmax over candidate rewards. The
probability measure $Q_{k}(\zeta_{j}^{k},\varrho_{j}^{k},d\zeta)$ is the
transition probability given candidate $\zeta_{j}^{k}$ and guidance signal
$\varrho_{j}^{k}$. Hence, $p_{k+1}(d\zeta)$ is a mixture of distributions
weighted by observed rewards in the current iteration $k$, which can be
sampled by the standard superposition method: at first the index $j$ is
sampled from the discrete distribution $\\{\nu_{j}^{k}\\}$, followed by
sampling from $Q_{k}(\zeta_{j}^{k},\varrho_{j}^{k},d\zeta)$. For example, in
our algorithm for CityLearn,

$Q_{k}(z,\varrho,d\zeta)\sim\exp(\|\zeta-z-\alpha_{k}\varrho\|/\iota_{k})\mu(d\zeta),$
(7)

where $\|\cdot\|$ is the Euclidean norm, $\mu(d\zeta)$ is a uniform measure
over $\mathcal{Z}$, and $\iota_{k}>0$ and $\alpha_{k}\geq 0$ are such that
their sum over time is bounded. Other candidates are possible and can still
ensure convergence to global optimal, as long as certain conditions are met;
intuitively, we require that the span of $Q_{k}$ decreases over time but not
so rapidly that it fails to reach a global optimum. Note that to simplify the
presentation, in the above algorithm, we assume that each candidate policy is
evaluated only on one episode; extending this to the case of multiple episodes
is straightforward (e.g., we would instead take the average of the evaluations
among the episodes in the computation of weights (6)).

##  4 Theoretical analysis

We now analyze the convergence property of the sequence generated by Algorithm
1. The following notations are used: $f(\zeta)=\mathbb{E}[R(\pi_{\zeta})]$ is
the expected episodic reward of policy $\pi_{\zeta}$,
$\Upsilon=\operatorname*{arg\,max}_{\zeta\in\mathcal{Z}}f(\zeta)$ is the set
of global maximizers (may not be unique),
$f^{*}=\max_{\zeta\in\mathcal{Z}}f(\zeta)$ is the global maximum, and
$\lambda(d\zeta)$ is some measure over $\Upsilon$;
$\mathbb{B}(\zeta,\epsilon)=\\{\zeta^{\prime}\in\mathcal{Z}:\|\zeta^{\prime}-\zeta\|\leq\epsilon\\}$
is a ball centered at $\zeta$ with radius $\epsilon$,
$\mathbb{B}^{*}(\epsilon)=\\{\zeta\in\mathcal{Z}:\min_{\zeta^{\prime}\in\Upsilon}\|\zeta^{\prime}-\zeta\|\leq\epsilon\\}$
is a set of points that are $\epsilon$ away from the optimal solution set
$\Upsilon$; $\delta_{\zeta}(dz)$ is the probability measure concentrated at
the point $\zeta$. We use $\Rightarrow$ to denote the weak convergence of
measures. We can consider $\|\cdot\|$ as any norm (e.g., Euclidean norm).

The measures $p_{k+1}(d\zeta)$, $k\in\mathbb{N}$ defined in (5) are
distributions of random points $\zeta_{j}^{k+1}$, for any $j\in[N_{k+1}]$,
conditional on the results of preceding evaluations of
$\\{R_{j}^{k}\\}_{j\in[N_{k}]}$ and realizations of
$\\{\zeta_{j}^{k},\varrho_{j}^{k},\xi_{j}^{k}\\}_{j\in[N_{k}]}$. Let
$P_{k}(d\zeta_{1},...,d\zeta_{N_{k}})$ represent their unconditional joint
distributions at iteration $k$, and

$\tilde{P}_{k}(d\zeta)=\int_{\mathcal{Z}^{N_{k}-1}}P_{k}(d\zeta,dz_{2},...,dz_{N_{k}})$

is the unconditional marginal distribution (note that we introduce $z$ for
$\zeta$ as the need arises in integration).

The formalism of the guidance signal requires some basics of random process
and measure theory (details can be found in the appendix). Essentially, the
guidance signal $\varrho_{j}^{k}\in\mathcal{Z}^{\prime}$ is a random variable
(adapted to the $\sigma$-algebra generated by the trajectory within an
episode) with probability measure $M_{k}(\zeta_{j}^{k},d\varrho)$. Note that
$M_{k}(\zeta_{j}^{k},d\varrho)$ is dependent on $\zeta_{j}^{k}$ because the
stochastic process that generates the trajectory depends on policy
$\pi_{\zeta_{j}^{k}}$, but $M_{k}(\zeta_{j}^{k},d\varrho)$ is conditionally
independent of all other candidates
$\\{\zeta_{j^{\prime}}^{k}\\}_{j^{\prime}\neq j}$. For analysis, we make the
following assumptions.

######  Assumption 2.

The followings hold:

  1. (a)

$R_{j}^{k}=f({\zeta_{j}^{k}})+\xi_{j}^{k}$, where
$\xi_{j}^{k}\overset{iid}{\sim}F_{k}(d\xi)$ for any $k\in\mathbb{N}$ are
independent with distribution $F_{k}(d\xi)$ bounded on a finite interval
$[-c_{\xi},c_{\xi}]$ and $\mathbb{E}\exp(\xi_{j}^{k})=1$;

  2. (b)

$|f(\zeta)|\leq c_{f}$ for all $\zeta\in\mathcal{Z}$ and $\mathcal{Z}$ is
compact;

  3. (c)

there exists $\epsilon>0$ such that $f$ is continuous on
$\mathbb{B}^{*}(\epsilon)$;

  4. (d)

$Q_{k}(z,\varrho,d\zeta)=q_{k}(z,\varrho,\zeta)\mu(d\zeta)$, with
$\sup_{z,\varrho,\zeta\in\mathcal{Z}}q_{k}(z,\varrho,\zeta)\leq L_{k}<\infty$
for all $k\in\mathbb{N}$, where $\mu$ is a probability measure on
$\mathcal{Z}$ such that $\mu(\mathbb{B}^{*}(\epsilon))>0$ for any
$\epsilon>0$; for any $z\in\mathcal{Z}$, the sequence of probability measures
$Q_{k}(z,\varrho,d\zeta)$ weakly converges to $\delta_{z}(d\zeta)$;

  5. (e)

$\\{N_{k}\\}$ is a sequence of natural numbers $N_{k}\in\mathbb{N}$ such that
$N_{k}\rightarrow\infty$ for $k\rightarrow\infty$;

  6. (f)

$\sup_{z,\varrho}M_{k}(z,d\varrho)<\infty$ for all $k\in\mathbb{N}$;

  7. (g)

$\tilde{P}_{1}(\mathbb{B}(\zeta,\epsilon))>0$ for all $\epsilon>0$ and
$\zeta\in\mathcal{Z}$;

  8. (h)

for any $\epsilon>0$, there are $\delta>0$ and a natural $\bar{k}$ such that
$\tilde{P}_{k}(\mathbb{B}^{*}(\epsilon))\geq\delta$ for all $k\geq\bar{k}$.

Let us comment on the assumptions above. Condition $(a)$ requires that the
evaluation noise be independent and bounded; the expectation requirement can
be satisfied for truncated log-normal distributions (Thompson 1950). The iid
requirement can be relaxed to mixing processes at the cost of more complex
analysis (Doukhan 2012); the boundedness condition, on the other hand, seems
necessary to keep iterates in the vicinity of global maximum if they are
already there. Condition $(b)$ is non-restrictive for practical problems.
Condition $(c)$ is natural since $\pi_{\zeta}$ is continuous by Lemma 1, and
can be satisfied if the true reward functions $r_{t}$ are continuous for all
$t\in[T]$. Assumptions $(d),(e),(f),(g)$, and $(h)$ formulate necessary
requirements on the parameters of the algorithm. Intuitively, conditions
$(d)$, $(e)$, and $(f)$ stipulate that the search becomes more “focused” over
time in order to concentrate on the global optima; however, conditions $(g)$
and $(h)$ indicates that the decrease of span cannot be too fast in order not
to miss the global optima. Condition $(f)$ on the guidance function $\phi$ can
be met by proper smoothing if needed. Condition $(e)$ can be relaxed to
$N_{k}=N$ for some finite integer $N$ for all $k\in\mathbb{N}$, but the
convergence will only be towards the vicinity of $\Lambda$ due to the finite
sample effect (see Lemma 4 in the appendix, which states the rate to be on the
order $N^{-1/2}$). Unlike $(c),(d),(e),(f)$, and $(g)$, condition $(h)$ is not
constructive; hence, we provide some verifiable conditions sufficient for
$(h)$ to hold in Corollary 1 (also see Corollary 2 in the appendix). Next, we
analyze the update rule of Algorithm 1.

######  Lemma 2.

The probability distribution $P_{k+1}(d\zeta_{1},...,d\zeta_{N_{k+1}})$ can be
written in terms of the distribution $P_{k}(d\zeta_{1},...,d\zeta_{N_{k}})$
as:

$\displaystyle\int_{\Omega^{N_{k}}}\\!\\!\chi_{k}(d\omega_{N_{k}})\\!\prod_{j=1}^{N_{k+1}}\\!\left\\{\beta(\omega_{N_{k}})\sum_{i=1}^{N_{k}}\Lambda(z_{i},\varrho_{i},\xi_{i},d\zeta_{j})\right\\},$
(8)

where $\Omega=\mathcal{Z}\times\mathcal{Z}^{\prime}\times[-c_{\xi},c_{\xi}]$,

$\displaystyle\omega_{N_{k}}=\\{z_{1},...,z_{N_{k}},\varrho_{1},...,\varrho_{N_{k}},\xi_{1},...,\xi_{N_{k}}\\}\in\Omega^{N_{k}},$
$\displaystyle\chi_{k}(d\omega_{N_{k}})=P_{k}(dz_{1},...,dz_{N_{k}})\prod_{j=1}^{N_{k}}F_{k}(d\xi_{j})M_{k}(z_{j},d\varrho_{j}),$
$\displaystyle\beta(\omega_{N_{k}})=\frac{1}{\sum_{j=1}^{N_{k}}\exp(f(z_{j})+\xi_{j})},\quad\text{
and}$
$\displaystyle\Lambda(z,\varrho,\xi,d\zeta)=\exp(f(z)+\xi)Q_{k}(z,\varrho,d\zeta).$

The proof is immediate by recognizing that the bracket term in (8) is the
conditional distribution $p_{k+1}(d\zeta_{j})$ defined in (5) and the
integration is over the distribution from the preceding iteration. We take the
product over $N_{k+1}$ candidates since they are drawn iid from $p_{k+1}$.

Now, we provide the main result on the convergence of $\tilde{P}_{k}(d\zeta)$
to some distribution $\lambda(d\zeta)$ over the global optima.

######  Theorem 1.

Suppose that Assumption 2 holds true, and let $\\{\tilde{P}_{k}\\}$ be the
sequence of unconditional marginal distributions determined by Algorithm 1.
Then, the distribution sequence weakly converges to some measure $\lambda$
over the optimal set, i.e., $\tilde{P}_{k}\Rightarrow\lambda$ as
$k\rightarrow\infty$.

The key stage of proof is to show that there exists a subsequence in
$\\{\tilde{P}_{k}\\}$ that weakly converges to

$\vartheta_{m}(d\zeta)=\frac{\exp(mf(\zeta))\kappa(d\zeta)}{\int\exp(mf(z))\kappa(dz)}$

for some measure $\kappa$, where $m$ is the index of the subsequence. The
above distribution is effectively a softmax function over function values and
converges to the extrema as $m\rightarrow\infty$.

All the conditions in Assumption 2 are natural with the exception of $(h)$,
which requires some further justification. In the following, we present a
sufficient condition for $(h)$ with a proper design of
$Q_{k}(z,\varrho,d\zeta)$, which applies to the case of noisy function
evaluations; see the appendix for another example in the case of noiseless
function evaluations.

######  Corollary 1.

Under Assumption 2 (except $(h)$), and let the transition probability
$Q_{k}(z,\varrho,d\zeta)$ be

$Q_{k}(z,\varrho,d\zeta)=c_{k}(z,\varrho)\psi((\zeta-z-\alpha_{k}\varrho)/\iota_{k})\mu(d\zeta),$
(9)

where
$c_{k}(z,\varrho)=({\int\psi((\zeta-z-\alpha_{k}\varrho)/\iota_{k})\mu(d\zeta)})^{-1}$
is the normalization term, $\psi$ is a continuous symmetrical finite density,
and

$\iota_{k}>0,\quad\sum_{k=1}^{\infty}\iota_{k}<\infty,\quad\alpha_{k}\geq
0,\quad\sum_{k=1}^{\infty}\alpha_{k}<\infty.$

Then, there exists a sequence of natural numbers $\\{N_{k}\\}$ such that
$\\{\tilde{P}_{k}\\}$ weakly converges to $\lambda$.

Our analysis accounts for the effect of trajectory-based guidance, which is a
novel contribution to the ES literature. By examining the proof, we can relax
the condition that $\sum_{k=1}^{\infty}\alpha_{k}<\infty$, i.e., continue
applying the guidance without the need to diminish its impact in the long run,
as long as the guidance signal “approximately” points to the global optima in
the proximity (see the appendix for exact conditions). However, such guidance
can be difficult to design or even verify in practice; thus, it is still
advisable to relinquish human knowledge and let data drive the decision,
eventually.

##  5 Results from the CityLearn Challenge

_Challenge overview._ The competition has an online setup with a simulation
period of 1 or 4 years, where agents exploit the best policies to optimize the
coordination strategy. The goal of each agent is to minimize environmental
costs, such as ramping costs, peak demands, 1-load factor, and carbon
emissions. The state space contains information such as daylight hours,
outdoor temperature, storage device state of charge (SOC), net electricity
consumption of the building, carbon intensity of the power grid, among a total
of 30 continuous states. The agent is allowed to control the
charging/discharging actions of storage devices for domestic hot water (DHW),
chilled water, and electricity (i.e., 3 continuous actions per building). The
environment is seen as a blackbox to the agent as a standard RL setup, where
the transition dynamics depend on the responses of various devices (e.g., air-
to-water heat pumps, electric heaters) as well as the energy loads of
buildings, which include space cooling, dehumidification, DHW demand, and
solar generation.

_Evaluation._ The submission of each team is evaluated on a set of metrics,
including: _(1)_ ramping: $\sum|e_{t}-e_{t-1}|$, where $e$ is the net
electricity consumption at each time step; _(2)_ 1-load factor: average net
electricity load divided by maximum electricity load; _(3)_ average daily peak
demand; _(4)_ maximum peak electricity demand; _(5)_ total electricity
consumed; _(6)_ carbon emissions. The competition evaluates performance by
computing the ratio of costs with respect to a rule-based controller
(RBC)—lower ratios indicate better performances.222Note that the RBC
controller is ubiquitous in traditional building control systems and is a
simple form of “take action $a_{h}$ in hour $h$,” where $a_{h}$ is a constant
independent of current states except for the hour of the day ($h\in[24]$).
The average of the above metrics for the full simulated period is the _total
score_ , while the average of the metrics _(1)_ -_(4)_ is _coordination
score_. The performance of the top 4 teams is listed in Table 1. Refer to
(Vazquez-Canteli et al. 2020) for more details on the contest.

Figure 2: Learning curves of ZO-iRL and baselines. We perform 10 runs on each
baseline to obtain performance plots with standard deviations for Climate
Zones 1 to 5.

ZO-iRL(ours) ICD-CA IDLab Breakfast Club Total score 0.944 1.070 1.070 1.130
Coord. score 0.915 1.107 1.098 1.095

Table 1: Total and coordination scores of the top 4 teams in 2021 CityLearn
Challenge.

SAC A2C DDPG DQN PPO TD3 MARLISA ZO-iRL ramping cost 1.145 (0.015) 1.189
(0.002) 1.174 (0.029) 1.302 (0.008) 1.638 (0.005) 1.178 (0.026) 1.022 (0.010)
0.781 (0.005) 1-load factor 1.158 (0.002) 1.146 (0.002) 1.143 (0.007) 1.159
(0.003) 1.168 (0.003) 1.142 (0.004) 1.026 (0.006) 1.010 (0.008) avg. daily
peak 1.180 (0.002) 1.184 (0.008) 1.195 (0.011) 1.212 (0.002) 1.242 (0.001)
1.193 (0.006) 1.015 (0.001) 0.996 (0.001) peak demand 1.077 (0.008) 1.088
(0.007) 1.100 (0.010) 1.115 (0.009) 1.132 (0.006) 1.098 (0.013) 1.000 (6e-5)
0.962 (0.005) net electric. peak 0.995 (0.001) 0.994 (0.002) 0.997 (8e-4)
0.997 (1e-4) 1.003 (7e-5) 0.997 (7e-4) 1.000 (5e-4) 1.006 (2e-4) carbon
emissions 1.000 (0.001) 1.000 (0.002) 1.003 (7e-4) 1.005 (1e-4) 1.009 (7e-5)
1.004 (7e-4) 1.001 (5e-4) 1.007 (2e-4) total score 1.092 (0.003) 1.101 (0.003)
1.102 (0.008) 1.132 (0.002) 1.199 (0.001) 1.102 (0.006) 1.011 (0.002) 0.962
(0.001)

Table 2: Comparison with baselines: SAC (Kathirgamanathan et al. 2020) and
MARLISA (Vazquez-Canteli, Henze, and Nagy 2020) have been officially
implemented for CityLearn; other baselines are implemented by (Raffin et al.
2019). The reported values are the average and standard deviation (in
brackets) across 10 independent runs on Climate Zone 1 data.

ZO-iRL: zeroth-order implicit RL. 333We name our method ZO-iRL because the
policy is implicitly determined by solving an optimization problem and the
learning algorithm is zeroth-order in an RL setting. As our method is designed
for single-agent episodic RL, we first reduce the original task that consists
of a single period of 1 or 4 years into episodes of 24 hours. We use the per-
step reward $-\max(0,e_{t})^{3}$ as recommended by (Vazquez-Canteli et al.
2020), where $e_{t}$ is the net electricity consumption (or generation if
$e_{t}<0$). This reward favours consumption patterns that are smoothly
averaged without demand peaks, aligned with multiple metrics used in the
evaluation, such as the 1-load factor and peak electricity demand. Another
reduction is from multi-agent RL to single-agent RL, where each building’s
policy is updated independently, reducing the problem to decentralized control
with additive rewards; such a reduction is computationally efficient for
large-scale problems (De Nijs et al. 2021). We omit the notational dependence
on candidate $j$ and iteration $k$ when presenting the method.

_Optimization planner._ We instantiate the optimization in (2) as follows. The
planned states $\bar{s}_{t}$ consist of state variables such as net
electricity consumption and SOCs of storage devices; the action
$\bar{a}_{t}\in\mathcal{A}$ is the action of the MDP; the surrogate reward

$\bar{r}_{t}(\bar{s}_{t};\zeta)=-|e_{t}-e_{t-1}|-\theta_{t}e_{t}$

is a combination of the negated ramping cost and the “virtual” electricity
cost, where $\zeta=\\{\theta_{t}\in[0,5]\\}_{t\in[24]}$ can be viewed as
virtual electricity prices to be learned to encourage desirable consumption
patterns (e.g., load flattening and smoothing). Intuitively, a higher value of
$\theta_{t}$ discourages planned electricity consumption in the corresponding
hour $t$.

The inequalities are grouped into technological constraints (e.g.,
maximum/minimum cooling power) and constraints on states and actions. The
equalities are grouped into physics accounting for energy balances (i.e.,
consumption equal to supply) and technology (e.g., SOC update rules). Further
details are provided in the appendix. Note that to set up the optimization
(2), we also need to provide predictions of energy demands and solar
generation. For simplicity, our predictors are based on a simple averaging
scheme that takes the average in the corresponding hours of the last 2 weeks
of data; thus, there are no specific parameters to learn.

_Transition and guidance._ We use (7) as the transition probability, with
variance $\iota_{k}=0.4/k^{2}$ that is initialized to 0.4 and decreases by
$k^{2}$ in each iteration. The guidance signal $\varrho$ is computed as
follows. By the end of each episode, we examine the net electricity usage in
the past 24 hours, $e_{t}$ for $t\in[24]$ and find the top 2 hours with the
most electricity usage, denoted by $t_{1}$ and $t_{2}$. Then, the guidance
signal $\varrho_{t}$ is $0.02$ if $t\in\\{t_{1},t_{2}\\}$ and $-0.04/22$
otherwise. Note that we have centered the signal ($\sum_{t}\varrho_{t}=0$) by
assigning negative values for hours other than peaks. We choose $\alpha_{k}=1$
for all $k$ over the entire 4-year period, as there is no training phase in
the CityLearn Challenge and we prefer to adapt quickly during the test phase;
this is not a violation of our theory, as we can choose to diminish
$\alpha_{k}$ after a while to still satisfy the condition
$\sum_{k=1}^{\infty}\alpha_{k}<\infty$.

Results. For baselines, we use the implementation of (Raffin et al. 2019) with
the default ADAM optimizer, where the policy is an NN architecture with tanh
activation and two layers of 256 units each. From Table 2, we see that ZO-iRL
has achieved the lowest cost ratios (i.e., best scores) of all, which is
consistent with the official result of the competition (Table 1). In
particular, as shown in Fig. 2, ZO-iRL is able to find a good policy in the
first few months, while baselines seem to struggle; we speculate that more
samples would eventually improve the performance of baselines, and all methods
may benefit from schemes to handle the potentially nonstationary environment
due to seasonal patterns.

##  6 Conclusion and future directions

We presented a novel adaptive optimization framework that has been shown to be
very effective for energy storage management. Using solution functions as
policies offers a promising way to introduce data-driven algorithms into the
real world where convex optimization has been widely adopted. To adapt the
optimization parameters, we developed an evolutionary search algorithm that
can incorporate insights from control trajectory data as guidance for
parameter updates. The method outperforms several baselines and ranked first
in the latest 2021 CityLearn Challenge. Some potential future directions could
be to extend the proposed framework to other methods such as Bayesian
optimization or first-order methods such as actor-critic.

##  7 Acknowledgments

We would like to thank all the ZORL/ROLEVT team members at the CityLearn
Challenge 2021: Qasim Wani, Zhiyao Zhang, and Mingyu Kim. We would also like
to thank the organizers of the CityLearn Challenge: Dr. José R. Vázquez-
Canteli, Dr. Zoltan Nagy, Dr. Gregor Henze, and Sourav Dey for organizing this
competition and inspiring novel RL-based approaches to the contemporary
challenges facing the power and energy sector. This project is funded by NSF,
the Commonwealth Cyber Initiative (CCI), C3.ai Digital Transformation
Institute, and the U.S. Department of Energy.

## References

  * Abedi, Gaudard, and Romerio (2019) Abedi, A.; Gaudard, L.; and Romerio, F. 2019.  Review of major approaches to analyze vulnerability in power system.  _Reliability Engineering & System Safety_, 183: 153–172.
  * Agrawal et al. (2019) Agrawal, A.; Amos, B.; Barratt, S.; Boyd, S.; Diamond, S.; and Kolter, J. Z. 2019\.  Differentiable convex optimization layers.  _Advances in neural information processing systems_ , 32.
  * Agrawal et al. (2020) Agrawal, A.; Barratt, S.; Boyd, S.; and Stellato, B. 2020.  Learning convex optimization control policies.  In _Learning for Dynamics and Control_ , 361–373. PMLR.
  * Berge (1997) Berge, C. 1997.  _Topological Spaces: including a treatment of multi-valued functions, vector spaces, and convexity_.  Courier Corporation.
  * Bertsekas (2019) Bertsekas, D. 2019.  _Reinforcement learning and optimal control_.  Athena Scientific.
  * Billingsley (2013) Billingsley, P. 2013.  _Convergence of probability measures_.  John Wiley & Sons.
  * Blum et al. (1958) Blum, J.; Chernoff, H.; Rosenblatt, M.; and Teicher, H. 1958.  Central limit theorems for interchangeable processes.  _Canadian Journal of Mathematics_ , 10: 222–229.
  * Borrelli, Bemporad, and Morari (2017) Borrelli, F.; Bemporad, A.; and Morari, M. 2017.  _Predictive control for linear and hybrid systems_.  Cambridge University Press.
  * Boyd, Boyd, and Vandenberghe (2004) Boyd, S.; Boyd, S. P.; and Vandenberghe, L. 2004.  _Convex optimization_.  Cambridge university press.
  * Chen et al. (2022) Chen, X.; Qu, G.; Tang, Y.; Low, S.; and Li, N. 2022.  Reinforcement learning for selective key applications in power systems: Recent advances and future challenges.  _IEEE Transactions on Smart Grid_.
  * Conti et al. (2018) Conti, E.; Madhavan, V.; Petroski Such, F.; Lehman, J.; Stanley, K.; and Clune, J. 2018.  Improving exploration in evolution strategies for deep reinforcement learning via a population of novelty-seeking agents.  _Advances in neural information processing systems_ , 31.
  * De Nijs et al. (2021) De Nijs, F.; Walraven, E.; De Weerdt, M.; and Spaan, M. 2021.  Constrained multiagent Markov decision processes: A taxonomy of problems and algorithms.  _Journal of Artificial Intelligence Research_ , 70: 955–1001.
  * Dempe and Zemkoho (2020) Dempe, S.; and Zemkoho, A. 2020.  _Bilevel optimization_.  Springer.
  * Dontchev and Rockafellar (2009) Dontchev, A. L.; and Rockafellar, R. T. 2009.  _Implicit functions and solution mappings_ , volume 543.  Springer.
  * Doukhan (2012) Doukhan, P. 2012.  _Mixing: properties and examples_ , volume 85.  Springer Science & Business Media.
  * Dulac-Arnold et al. (2021) Dulac-Arnold, G.; Levine, N.; Mankowitz, D. J.; Li, J.; Paduraru, C.; Gowal, S.; and Hester, T. 2021.  Challenges of real-world reinforcement learning: definitions, benchmarks and analysis.  _Machine Learning_ , 110(9): 2419–2468.
  * Facchinei and Pang (2007) Facchinei, F.; and Pang, J.-S. 2007.  _Finite-dimensional variational inequalities and complementarity problems_.  Springer Science & Business Media.
  * Frazier (2018) Frazier, P. I. 2018.  Bayesian optimization.  In _Recent Advances in Optimization and Modeling of Contemporary Problems_ , 255–278. INFORMS.
  * Gangwani and Peng (2018) Gangwani, T.; and Peng, J. 2018.  Policy Optimization by Genetic Distillation.  In _International Conference on Learning Representations_.
  * Ghadimi and Lan (2013) Ghadimi, S.; and Lan, G. 2013.  Stochastic first-and zeroth-order methods for nonconvex stochastic programming.  _SIAM Journal on Optimization_ , 23(4): 2341–2368.
  * Ghadimi, Perkins, and Powell (2020) Ghadimi, S.; Perkins, R. T.; and Powell, W. B. 2020.  Reinforcement Learning via Parametric Cost Function Approximation for Multistage Stochastic Programming.  _arXiv preprint arXiv:2001.00831_.
  * Guez et al. (2018) Guez, A.; Weber, T.; Antonoglou, I.; Simonyan, K.; Vinyals, O.; Wierstra, D.; Munos, R.; and Silver, D. 2018.  Learning to search with mctsnets.  In _International conference on machine learning_ , 1822–1831. PMLR.
  * Hajek (2015) Hajek, B. 2015.  _Random processes for engineers_.  Cambridge university press.
  * Haydari and Yilmaz (2020) Haydari, A.; and Yilmaz, Y. 2020.  Deep reinforcement learning for intelligent transportation systems: A survey.  _IEEE Transactions on Intelligent Transportation Systems_.
  * Jin et al. (2023) Jin, M.; Khattar, V.; Kaushik, H.; Sel, B.; and Jia, R. 2023.  On Solution Functions of Optimization: Universal Approximation and Covering Number Bounds.  In _AAAI Conference on Artificial Intelligence (AAAI)_.
  * Karkus, Hsu, and Lee (2017) Karkus, P.; Hsu, D.; and Lee, W. S. 2017.  Qmdp-net: Deep learning for planning under partial observability.  _Advances in neural information processing systems_ , 30.
  * Kathirgamanathan et al. (2020) Kathirgamanathan, A.; Twardowski, K.; Mangina, E.; and Finn, D. P. 2020.  A Centralised Soft Actor Critic Deep Reinforcement Learning Approach to District Demand Side Management through CityLearn.  In _Proceedings of the 1st International Workshop on Reinforcement Learning for Energy Management in Buildings & Cities_, 11–14.
  * Khadka and Tumer (2018) Khadka, S.; and Tumer, K. 2018.  Evolution-guided policy gradient in reinforcement learning.  _Advances in Neural Information Processing Systems_ , 31.
  * Liu et al. (2020) Liu, S.; Chen, P.-Y.; Kailkhura, B.; Zhang, G.; Hero III, A. O.; and Varshney, P. K. 2020.  A primer on zeroth-order optimization in signal processing and machine learning: Principals, recent advances, and applications.  _IEEE Signal Processing Magazine_ , 37(5): 43–54.
  * Mania, Guy, and Recht (2018) Mania, H.; Guy, A.; and Recht, B. 2018.  Simple random search of static linear policies is competitive for reinforcement learning.  _Advances in Neural Information Processing Systems_ , 31.
  * Moerland, Broekens, and Jonker (2020) Moerland, T. M.; Broekens, J.; and Jonker, C. M. 2020.  Model-based reinforcement learning: A survey.  _arXiv preprint arXiv:2006.16712_.
  * Nagy et al. (2021) Nagy, Z.; Vázquez-Canteli, J. R.; Dey, S.; and Henze, G. 2021.  The CityLearn Challenge 2021.  In _Proceedings of the 8th ACM International Conference on Systems for Energy-Efficient Buildings, Cities, and Transportation_ , 218–219.
  * Nian, Liu, and Huang (2020) Nian, R.; Liu, J.; and Huang, B. 2020.  A review on reinforcement learning: Introduction and applications in industrial process control.  _Computers & Chemical Engineering_, 139: 106886.
  * Pflug and Pichler (2014) Pflug, G. C.; and Pichler, A. 2014.  _Multistage stochastic optimization_ , volume 1104.  Springer.
  * Powell (2020) Powell, W. 2020.  Reinforcement Learning and Stochastic Optimization: A Unified Framework for Sequential Decisions.  _Princeton NJ_.
  * Prakash et al. (2020) Prakash, B.; Waytowich, N.; Ganesan, A.; Oates, T.; and Mohsenin, T. 2020.  Guiding safe reinforcement learning policies using structured language constraints.  _UMBC Student Collection_.
  * Prékopa (2013) Prékopa, A. 2013.  _Stochastic programming_ , volume 324.  Springer Science & Business Media.
  * Racanière et al. (2017) Racanière, S.; Weber, T.; Reichert, D.; Buesing, L.; Guez, A.; Jimenez Rezende, D.; Puigdomènech Badia, A.; Vinyals, O.; Heess, N.; Li, Y.; et al. 2017.  Imagination-augmented agents for deep reinforcement learning.  _Advances in neural information processing systems_ , 30.
  * Raffin et al. (2019) Raffin, A.; Hill, A.; Ernestus, M.; Gleave, A.; Kanervisto, A.; and Dormann, N. 2019\.  Stable baselines3.
  * Rockafellar and Wets (2009) Rockafellar, R. T.; and Wets, R. J.-B. 2009.  _Variational analysis_ , volume 317.  Springer Science & Business Media.
  * Rolnick et al. (2022) Rolnick, D.; Donti, P. L.; Kaack, L. H.; Kochanski, K.; Lacoste, A.; Sankaran, K.; Ross, A. S.; Milojevic-Dupont, N.; Jaques, N.; Waldman-Brown, A.; et al. 2022\.  Tackling climate change with machine learning.  _ACM Computing Surveys (CSUR)_ , 55(2): 1–96.
  * Salimans et al. (2017) Salimans, T.; Ho, J.; Chen, X.; Sidor, S.; and Sutskever, I. 2017.  Evolution strategies as a scalable alternative to reinforcement learning.  _arXiv preprint arXiv:1703.03864_.
  * Schrittwieser et al. (2020) Schrittwieser, J.; Antonoglou, I.; Hubert, T.; Simonyan, K.; Sifre, L.; Schmitt, S.; Guez, A.; Lockhart, E.; Hassabis, D.; Graepel, T.; et al. 2020.  Mastering atari, go, chess and shogi by planning with a learned model.  _Nature_ , 588(7839): 604–609.
  * Snoek, Larochelle, and Adams (2012) Snoek, J.; Larochelle, H.; and Adams, R. P. 2012.  Practical bayesian optimization of machine learning algorithms.  _Advances in neural information processing systems_ , 25.
  * Spall (2005) Spall, J. C. 2005.  _Introduction to stochastic search and optimization: estimation, simulation, and control_ , volume 65.  John Wiley & Sons.
  * Szita and Lörincz (2006) Szita, I.; and Lörincz, A. 2006.  Learning Tetris using the noisy cross-entropy method.  _Neural computation_ , 18(12): 2936–2941.
  * Tamar et al. (2016) Tamar, A.; Wu, Y.; Thomas, G.; Levine, S.; and Abbeel, P. 2016.  Value iteration networks.  _Advances in neural information processing systems_ , 29.
  * Thompson (1950) Thompson, H. 1950.  Truncated normal distributions.  _Nature_ , 165(4194): 444–445.
  * Vazquez-Canteli et al. (2020) Vazquez-Canteli, J. R.; Dey, S.; Henze, G.; and Nagy, Z. 2020.  CityLearn: Standardizing Research in Multi-Agent Reinforcement Learning for Demand Response and Urban Energy Management.  _arXiv preprint arXiv:2012.10504_.
  * Vazquez-Canteli, Henze, and Nagy (2020) Vazquez-Canteli, J. R.; Henze, G.; and Nagy, Z. 2020.  MARLISA: Multi-Agent Reinforcement Learning with Iterative Sequential Action Selection for Load Shaping of Grid-Interactive Connected Buildings.  In _Proceedings of the 7th ACM International Conference on Systems for Energy-Efficient Buildings, Cities, and Transportation_ , 170–179.
  * Zhigljavsky (2012) Zhigljavsky, A. A. 2012.  _Theory of global random search_ , volume 65.  Springer Science & Business Media.
  * Zhou, Yu, and Qian (2019) Zhou, Z.-H.; Yu, Y.; and Qian, C. 2019.  _Evolutionary learning: Advances in theories and algorithms_.  Springer.

##  Appendix A Proof of results in the main paper

###  A.1 Formalism of the guidance signal

The formalism of the guidance signal requires some basics from random process
and measure theory (interested readers are referred to (Hajek 2015)). We keep
our presentation minimal but sufficient enough to carry out the analysis.
Consider a stochastic process $(S_{t}^{j,k},A_{t}^{j,k})_{t\in[T]}$ defined by
policy $\pi_{\zeta_{j}^{k}}$ interacting with the MDP environment, where
$S_{t}^{j,k}$ and $A_{t}^{j,k}$ are random variables representing the state
and action at time $t$. Let
$\mathcal{I}_{t}=(\mathcal{S}\times\mathcal{A})^{t}$ be the set of possible
histories up to time step $t$ within an episode, and
$I_{t}^{j,k}\coloneqq(S_{1}^{j,k},A_{1}^{j,k},...,S_{t}^{j,k},A_{t}^{j,k})\in\mathcal{I}_{t}$
is a random vector taking values in $\mathcal{I}_{t}$ containing all state-
action pairs observed up to step $t$. Denote by $\mathcal{F}_{t}^{j,k}$ a non-
decreasing sequence of $\sigma$-algebra (a _filtration_) generated by
$I_{t}^{j,k}$. Then, the guidance signal
$\varrho_{j}^{k}\in\mathcal{Z}^{\prime}$ is a random variable adapted to the
filtration $\mathcal{F}_{T}^{j,k}$, i.e., $\varrho_{j}^{k}$ is
$\mathcal{F}_{T}^{j,k}$-measurable, with probability measure
$M_{k}(\zeta_{j}^{k},d\varrho)$ associated with a properly defined probability
space, the existence of which is ensured by the Ionescu-Tulcea theorem. Note
that the distribution of $\varrho_{j}^{k}$ depends on $\zeta_{j}^{k}$, since
the stochastic process $I_{t}^{j,k}$ is determined by the policy
$\pi_{\zeta_{j}^{k}}$, but is conditionally independent of all other
candidates $\zeta_{j^{\prime}}^{k}$ for $j^{\prime}\neq j$.

###  A.2 Proof of Lemma 1

Let $\Phi(s_{t},\zeta)$ represent the feasible set of (2). By Assumption 1,
$\Phi(s_{t},\zeta)$ is convex for fixed $s_{t}$ and $\zeta$ and has a nonempty
interior. This implies that $\Phi(s_{t},\zeta)$ is continuous in $s_{t}$ and
$\zeta$ (Rockafellar and Wets 2009, example 5.10). Hence, by Berge maximum
theorem (Berge 1997), $\pi_{\zeta}(s_{t})$ is upper hemicontinuous in $\zeta$
for fixed $s_{t}\in\mathcal{S}$. However, we know that $\pi_{\zeta}(s_{t})$
contains a single point due to the strict convexity of the objective function.
Thus, for fixed $s_{t}\in\mathcal{S}$, $\pi_{\zeta}(s_{t})$ is a single-valued
function continuous in its parameter $\zeta$.

###  A.3 Proof of Theorem 1

Select from $\\{\tilde{P}_{k}\\}$ a weakly convergent subsequence
$\\{\tilde{P}_{k_{i}}\\}$, which is possible due to Prohorov’s theorem
(Billingsley 2013, Ch. 6), and denote the limit by $\kappa(d\zeta)$. By Lemma
4, we have that

$\tilde{P}_{k+1}(d\zeta)=\left(\int\tilde{P}_{k}(dz)\exp(f(z))\right)^{-1}\int\tilde{P}_{k}(dz)\exp(f(z))\Big{(}Q_{k}(z,\varrho,d\zeta)M_{k}(z,d\varrho)+\Delta_{N_{k}}(d\zeta)\Big{)}.$
(10)

By Assumption 2 $(d)$ and $(h)$, it follows that the subsequence
$\\{\tilde{P}_{k_{i}+1}\\}$ weakly converges to the distribution
$\vartheta_{1}(d\zeta)=c_{1}\exp(f(\zeta))\kappa(d\zeta)$, where $c_{1}$ is
the normalization constant. Similarly, the subsequence
$\\{\tilde{P}_{k_{i}+m}\\}$ weakly converges to the distribution

$\vartheta_{m}(d\zeta)=\frac{\exp(mf(\zeta))\kappa(d\zeta)}{\int\exp(mf(z))\kappa(dz)},$

which, by Lemma 3, converges to $\lambda$. Thus, by the standard
diagonalization argument (Billingsley 2013), we can show that there exists a
subsequence $\\{\tilde{P}_{k_{j}}\\}$ that weakly converges to $\lambda$.
Applying Lemma 4 again yields that $\\{\tilde{P}_{k_{j}+1}\\}$ converges to
the same limit. Thus, any subsequence of $\\{\tilde{P}_{k}\\}$ converges to
this limit, and the same holds for the sequence itself.

###  A.4 Proof of Corollary 1

Under Assumption 2 (except $(h)$), the distributions (16) have continuous
densities with respect to the Lebesgue measure. Let
$A(\epsilon)=\\{\zeta\in\mathcal{Z}:f(\zeta)\geq f^{*}-\epsilon\\}$. By (9)
and Lemma 2, we have that $\tilde{P}_{k}(d\zeta)>0$ for any $k\in\mathbb{N}$.
Fix an arbitrary $\delta>0$. We shall choose $\\{N_{k}\\}$ such that for any
$k\geq k_{n}$ and $\epsilon>0$, the following holds

$\tilde{P}_{k+1}(A(\epsilon+\epsilon_{k}))\geq(1-\delta_{k})\tilde{P}_{k}(A(\epsilon)),$
(11)

where

$0<\delta_{k}<1\quad\text{for
$k\in\mathbb{N}$},\qquad\quad\sum_{k\in\mathbb{N}}\delta_{k}<\infty,$ (12)

and $\epsilon_{k}\geq 0$ are determined in terms of $\iota_{k}$, $\alpha_{k}$,
and the sizes of the support of density $\psi$,

$\sum_{k=1}^{\infty}\epsilon_{k}\leq\mathrm{constant}\sum_{k=1}^{\infty}\iota_{k}<\infty.$
(13)

Such sequence of $\\{N_{k}\\}$ and $k_{n}$ exist by Lemma 2, the finiteness of
$\psi$, and the condition that $\sum_{k=1}^{\infty}\alpha_{k}<\infty$. Next,
select $k_{o}\geq k_{n}$ such that

$\sum_{k=k_{o}}^{\infty}\epsilon_{k}<\frac{1}{2}\delta,$

and let $\delta_{1}=\tilde{P}_{k_{o}}(A(\delta/2)).$ Then, for any $k\geq
k_{o}$, we have

$\displaystyle\tilde{P}_{k+1}(A(\delta))$
$\displaystyle\geq\tilde{P}_{k_{o}}(A(\delta/2+\sum_{i=k_{o}}^{k}\delta_{i}))\prod_{i=k_{o}}^{k}(1-\delta_{i})$
$\displaystyle\geq\delta_{1}\prod_{i=k_{o}}^{\infty}(1-\delta_{i})$
$\displaystyle>0$

where the last inequality is implied by (12). The proof is complete.

Remarks on the guidance signal. From the proof of Corollary 1, it can be
observed that we can relax the condition that
$\sum_{k=1}^{\infty}\alpha_{k}<\infty$ as long as the guidance signal
$\alpha_{k}\varrho^{k}$ is chosen in such a way that (11), (12), and (13) are
satisfied. This means that we can continue applying the guidance signal
without the need to diminish its impact in the long run. However, (11) is
difficult to ensure, as it requires designing a guidance that always points to
the global optimal. Therefore, in practice, it is recommended to diminish the
effect of guidance and eventually let the data drive the decision.

###  A.5 An example of transition probability with noiseless function
evaluations

######  Corollary 2.

Under Assumption 2 (except for $(h)$), and further assume that $f$ can be
evaluated without noise (i.e., $\xi=0$). Let the transition probability
$Q_{k}(z,\varrho,A)$ be defined by

$\displaystyle Q_{k}(z,\varrho,A)=$ $\displaystyle\int 1_{\left\\{\zeta\in
A,f(z)\leq f(\zeta)\right\\}}T_{k}(z,\varrho,d\zeta)$
$\displaystyle+1_{\left\\{z\in A\right\\}}\int
1_{\left\\{f(\zeta)<f(z)\right\\}}T_{k}(z,\varrho,d\zeta),$ (14)

where $\\{T_{k}(z,\varrho,d\zeta)\\}$ weakly converges to $\delta_{z}(d\zeta)$
for all $z,\varrho\in\mathcal{Z}$. Then, there exists a sequence of natural
numbers $N_{k}$ such that the sequence of distributions $\\{\tilde{P}_{k}\\}$
weakly converges to $\lambda$ for $k\rightarrow\infty$.

###### Proof.

By Assumption 2 $(c)$ and $(g)$, we have that
$\tilde{P}_{1}(\mathbb{B}^{*}(\epsilon))>0$ for any $\epsilon>0$. By (2), we
have that

$\tilde{P}_{k}(\mathbb{B}^{*}(\epsilon))\geq\cdots\geq\tilde{P}_{1}(\mathbb{B}^{*}(\epsilon))>0$

for all $k\in\mathbb{N}$. Hence, Assumption 2 $(h)$ is satisfied. By Theorem
1, the claim is proved. ∎

Remarks. To implement the transition of (2), one first needs to sample a
variable $\zeta$ according to $T_{k}(z,\varrho,d\zeta)$ and observe its reward
value $f(\zeta)$; then, the output is $\zeta$ if $f(\zeta)\geq f(z)$ and $z$
otherwise. Such a scheme depends crucially on a reliable way of comparing
candidates (e.g., noiseless evaluation).

###  A.6 Supporting lemmas

######  Lemma 3.

Under Assumption 2 $(b),(c)$, and $(d)$, the sequence of distributions

$\frac{\exp(kf(\zeta))\mu(d\zeta)}{\int\exp(kf(z))\mu(dz)}\Rightarrow\lambda(d\zeta),$

i.e., weakly converges to $\lambda(d\zeta)$ for $k\rightarrow\infty$.

###### Proof.

By the definition of weak convergence, it suffices to show that for any
function $\Psi(\zeta)$ continuous on $\mathcal{Z}$, it holds that

$\lim_{k\rightarrow\infty}c_{k}\int\exp(kf(\zeta))\Psi(\zeta)\mu(d\zeta)=\int\Psi(\zeta)\lambda(d\zeta),$
(15)

where $c_{k}=1/\int\exp(kf(z))\mu(dz)$. To proceed, Let
$\mathbb{B}_{i}=\mathbb{B}(\epsilon_{i})=\\{\zeta\in\mathcal{Z}:\min_{\zeta^{\prime}\in\Lambda}\|\zeta^{\prime}-\zeta\|\leq\epsilon_{i}\\}$
and
$\mathbb{D}_{i}=\\{\zeta\in\mathcal{Z}:\min_{\zeta^{\prime}\in\Lambda}\|\zeta^{\prime}-\zeta\|\geq\epsilon_{i}\\}$,
for $i=0,1,2$ and some $\epsilon_{0},\epsilon_{1},$ and $\epsilon_{2}$ to be
determined. For any $\delta>0$, by continuity of $\Psi$, there exists
$\epsilon_{0}>0$ such that
$|\Psi(z)-\int\Psi(\zeta)\lambda(d\zeta)|\leq\delta$ for all
$z\in\mathbb{B}_{0}$. Choose some $\epsilon_{1}>0$ such that
$\epsilon_{1}<\epsilon_{0}$. Then, we have

$\displaystyle\left|c_{k}\int\exp(kf(\zeta))\Psi(\zeta)\mu(d\zeta)-\int\Psi(\zeta)\lambda(d\zeta)\right|$
$\displaystyle\leq
c_{k}\int_{\mathbb{B}_{1}}\exp(kf(z))\left|\Psi(z)-\int\Psi(\zeta)\lambda(d\zeta)\right|\mu(dz)+c_{k}\int_{\mathbb{D}_{1}}\exp(kf(z))\left|\Psi(z)-\int\Psi(\zeta)\lambda(d\zeta)\right|\mu(dz)$
$\displaystyle\leq\delta\underbrace{c_{k}\int_{\mathbb{B}_{1}}\exp(kf(z))\mu(dz)}_{(i)}+2\|\Psi\|_{\infty}\underbrace{c_{k}\int_{\mathbb{D}_{1}}\exp(kf(z))\mu(dz)}_{(ii)},$

where the first inequality is due to triangle inequality, and the second
inequality is due to the choice of $\epsilon_{1}$ (also, recall that
$\|\Psi\|_{\infty}=\sup|\Psi(z)|$). Hence, the lemma is proved if we can show
that $(i)\rightarrow 1$ and $(ii)\rightarrow 0$ as $k\rightarrow\infty$.

To this end, let $C_{1}=\sup_{\zeta\in\mathbb{D}_{1}}f(\zeta)$. By Assumption
2 $(c)$, there exists $\epsilon_{2}$ such that $0<\epsilon_{2}<\epsilon_{1}$,
and

$C_{2}=\inf_{\zeta\in\mathbb{B}_{2}}f(\zeta)>C_{1}.$

For any $k>0$, we have

$\int_{\mathbb{B}_{1}}\exp(kf(z)-kC_{1})\mu(dz)>\int_{\mathbb{B}_{2}}\exp(kf(z)-kC_{1})\mu(dz)\geq\int_{\mathbb{B}_{2}}\exp(k(C_{2}-C_{1}))\mu(dz).$

Thus,

$\frac{\int_{\mathbb{D}_{1}}\mu(dz)}{\int_{\mathbb{B}_{2}}\exp(k(C_{2}-C_{1}))\mu(dz)}\geq\underbrace{\frac{\int_{\mathbb{D}_{1}}\exp(kf(z))\mu(dz)}{\int_{\mathbb{B}_{1}}\exp(kf(z))\mu(dz)}}_{(iii)}\geq
0.$

By driving $k\rightarrow\infty$ to the limit and using the sandwich theorem,
we have that $(iii)\rightarrow 0$. This immediately implies that
$(i)\rightarrow 1$ and $(ii)\rightarrow 0$ as $k\rightarrow\infty$, hence
concluding the proof. ∎

######  Lemma 4.

Let Assumption 2 $(a),(b)$, and $(d)$ be fulfilled. Then, the marginal
distributions can be written as

$\displaystyle\tilde{P}_{k+1}(d\zeta)$
$\displaystyle=\left(\int\tilde{P}_{k}(dz)\exp(f(z))\right)^{-1}\int\tilde{P}_{k}(dz)\exp(f(z))Q_{k}(z,\varrho,d\zeta)M_{k}(z,d\varrho)+\Delta_{N_{k}}(d\zeta),$
(16)

where the signed measures $\Delta_{N_{k}}(d\zeta)$ converge to zero in
variation for ${N_{k}}\rightarrow\infty$ with the rate ${N_{k}}^{-1/2}.$

###### Proof.

For notational simplicity, we use $N$ for ${N_{k}}$ throughout the proof. By
Assumption 2 $(a)$ and Lemma 2, the marginal distribution
$\tilde{P}_{k+1}(d\zeta)$ is given by:

$\displaystyle\tilde{P}_{k+1}(d\zeta)$
$\displaystyle=\int_{\Omega^{N}}\chi_{k}(d\omega_{N})\left\\{\beta(\omega_{N})\sum_{i=1}^{N}\Lambda(\zeta_{i},\varrho_{i},\xi_{i},d\zeta)\right\\}$
$\displaystyle=\sum_{i=1}^{N}\int_{\Omega^{N}}\chi_{k}(d\omega_{N})\beta(\omega_{N})\Lambda(\zeta_{i},\varrho_{i},\xi_{i},d\zeta)$
$\displaystyle=\int_{\Omega^{N}}\chi_{k}(d\omega_{N})\left\\{N\beta(\omega_{N})\right\\}\Lambda(\zeta_{1},\varrho_{1},\xi_{1},d\zeta).$

which can be represented in the form of (16) with

$\displaystyle\Delta_{N}(d\zeta)$
$\displaystyle={\int_{\Omega^{N}}\chi_{k}(d\omega_{N})\Lambda(\zeta_{1},\varrho_{1},\xi_{1},d\zeta)\left\\{N\beta(\omega_{N})-\left(\int\tilde{P}_{k}(dz)\exp(f(z))\right)^{-1}\right\\}}$
$\displaystyle\quad+{\left(\int\tilde{P}_{k}(dz)\exp(f(z))\right)^{-1}\left\\{\int_{\Omega^{N}}\chi_{k}(d\omega_{N})\Lambda(\zeta_{1},\varrho_{1},\xi_{1},d\zeta)-\int_{\Omega}\tilde{P}_{k}(dz)\exp(f(z))Q_{k}(z,\varrho,d\zeta)M_{k}(z,d\varrho)\right\\}}$
$\displaystyle=(i)+(ii)$

We shall show that $(i)\rightarrow 0$ in variation for $N\rightarrow\infty$
and $(ii)=0$. Due to Assumption 2 $(d)$, the convergence of $(i)$ is
equivalent to the fact that $\int|v_{N}(\zeta)|\mu(d\zeta)\rightarrow 0$,
where

$v_{N}(z)={\int_{\Omega^{N}}\chi_{k}(d\omega_{N})\exp(f(\zeta_{1})+\xi_{1})q_{k}(\zeta_{1},\varrho_{1},z)\left\\{N\beta(\omega_{N})-\left(\int\tilde{P}_{k}(dz)\exp(f(z))\right)^{-1}\right\\}}.$

To proceed, let
$\gamma_{N}=\frac{1}{N}\sum_{i=1}^{N}\exp(f(\zeta_{i})+\xi_{i})$ and
$\psi(z)=\exp(f(\zeta_{1})+\xi_{1})q_{k}(\zeta_{1},\varrho_{1},z)$. Due to the
symmetrical dependence of random elements $\zeta_{1},...,\zeta_{N}$ and
$\varrho_{1},...,\varrho_{N}$, as well as the independence of
$\xi_{1},...,\xi_{N}$, the random variables $\gamma_{N}$ converge in mean for
$N\rightarrow\infty$ to some random variable $\gamma$ in dependent of all
$\gamma_{i}(\omega_{i})$, $y_{i}=f(\zeta_{i})+\xi_{i}$, for $i\in\mathbb{N}$,
and

$\mathbb{E}\gamma=\mathbb{E}\exp(y_{i})=\int\exp(f(\zeta)+\xi)\tilde{P}_{k}(d\zeta)F_{k}(d\xi).$

Equivalently, for any $\delta_{1}>0$, there exists $N_{\gamma}(\delta_{1})\geq
1$ such that $\mathbb{E}|\gamma_{N}-\gamma|<\delta_{1}$ for all $N\geq
N_{\gamma}(\delta_{1})$. Then,

$\displaystyle|v_{N}(z)|$
$\displaystyle=\left|\mathbb{E}\left(\frac{\psi(z)}{\gamma_{N}}\right)-\frac{\mathbb{E}\psi(z)}{\mathbb{E}\gamma}\right|$
(17)
$\displaystyle=\frac{1}{\mathbb{E}\gamma}\left|\mathbb{E}\left(\frac{\psi(z)\gamma}{\gamma_{N}}\right)-{\mathbb{E}\psi(z)}\right|$
(18)
$\displaystyle\leq\exp(c_{f})\left|\mathbb{E}\left(\frac{\psi(z)|\gamma-\gamma_{N}|}{\gamma_{N}}\right)\right|$
(19)
$\displaystyle\leq\exp(2c_{f})\|\psi\|_{\infty}\mathbb{E}|\gamma-\gamma_{N}|$
(20) $\displaystyle\leq
L_{k}\exp(3c_{f}+c_{\xi})\mathbb{E}|\gamma-\gamma_{N}|,$ (21)

where the second equality is due to the independence of $\gamma$ from
$\gamma_{N}$ and $\psi$, the first and second inequalities are due to
$\gamma,\gamma_{N}\geq\exp(-c_{f})$ (by Assumption 2 $(b)$), and the last
relation is due to $\|\psi\|_{\infty}\leq\exp(f(\zeta)+\xi))L_{k}\leq
L_{k}\exp(c_{f}+c_{\xi})$. In order to show that
$\int|v_{N}(z)|\mu(dz)\rightarrow 0$, we need to prove that for any $\delta>0$
and $z\in\mathcal{Z}$, there exists $N^{\star}(\delta,z)$ such that for $N\geq
N^{\star}(\delta,z)$, there holds $|v_{N}(z)|\leq\delta$. This can hold if one
takes $\delta_{1}=\delta L_{k}^{-1}\exp(-3c_{f}-c_{\xi})$ and
$N^{\star}(\delta,z)=N_{\gamma}(\delta_{1})$.

Now, by (21), we have that $\int|v_{N}(\zeta)|\mu(d\zeta)\leq
L_{k}\exp(3c_{f}+c_{\xi})\mathbb{E}|\gamma-\gamma_{N}|$. From the central
limit theorem for symmetrically dependent random variables (see (Blum et al.
1958)), it follows that $\mathbb{E}|\gamma-\gamma_{N}|=\mathcal{O}(N^{-1/2})$.
Consequently, we have shown that
$\int|v_{N}(\zeta)|\mu(d\zeta)=\mathcal{O}(N^{-1/2})$.

To show that $(ii)=0$, note that

$\displaystyle\int_{\Omega^{N}}\chi_{k}(d\omega_{N})\Lambda(\zeta_{1},\varrho_{1},\xi_{1},d\zeta)-\int_{\Omega}\tilde{P}_{k}(dz)\exp(f(z))Q_{k}(z,\varrho,d\zeta)M_{k}(z,d\varrho)$
$\displaystyle=\int_{\mathcal{Z}}\tilde{P}_{k}(dz)\exp(f(z))Q_{k}(z,\varrho,d\zeta)M_{k}(z,d\varrho)\left\\{\int\exp\left(\xi\right)F_{k}(d\xi)-1\right\\},$

which is $0$ by Assumption 2 $(a)$. Hence, we have concluded the proof. ∎

##  Appendix B Additional details for the CityLearn Challenge

###  B.1 Details of optimization model

We refer the reader to (Vazquez-Canteli et al. 2020) and the corresponding
online documentation444link: https://sites.google.com/view/citylearnchallenge
for the detailed setup of the contest. We will focus only on our strategy in
this document. In particular, we provide details on the construction of the
optimization model in 2. Denote the hourly index by $r\in\\{1,2,\cdots,T\\}$,
where $T=24$. Suppose we are at the beginning of the hour $r$. Then we need to
plan for the actions for the future hours until the end of the day and execute
the plan for the next hour $r$, a.k.a., rolling-horizon planning. Next, we
describe hyperparameters, variables, objective, and constraints in 2.

Hyperparameters. Hyperparameters are required to instantiate an optimization
and are not part of the optimization variables to be solved by an optimization
algorithm.

  * •

The hyperparameters to be set by prior knowledge include: _(1)_ electric
heater: efficiency $\eta_{\text{ehH}}$, nominal power
$E_{\text{max}}^{\text{ehH}}$; _(2)_ heat pump: technical efficiency
$\eta^{\text{hp}}_{\text{tech}}$, target cooling temperature
$t_{c}^{\text{hp}}$, nominal power $E_{\text{max}}^{\text{hpc}}$; _(3)_
electric battery: rate of decay $Cf^{\text{bat}}$, capacity $Cp^{\text{bat}}$,
efficiency $\eta_{t}^{\text{bat}}$; _(4)_ heat storage: rate of decay
$Cf^{\text{Hsto}}$, capacity $Cp^{\text{Hsto}}$, efficiency
$\eta_{t}^{\text{Hsto}}$; _(5)_ cooling storage: rate of decay
$Cf^{\text{Csto}}$, capacity $Cp^{\text{Csto}}$, efficiency
$\eta_{t}^{\text{Csto}}$.

  * •

The hyperparameters provided by the predictors include: _(1)_ hourly
coefficient of performance (COP) of heat pump
$\text{COP}^{C}_{t}=\eta_{\text{tech}}^{\text{hp}}\frac{t_{c}^{\text{hp}}+273.15}{\text{temp}_{t}-t_{c}^{\text{hp}}}$
, where $\text{temp}_{t}$ is the predicted outside temperature for hour $t$;
_(2)_ solar generation $E^{PV}_{t}$; _(3)_ electricity non-shiftable load
$E^{NS}_{t}$; _(4)_ heating demand $H^{bd}_{t}$; and _(5)_ cooling demand
$C^{bd}_{t}$. At hour $r$, the above predictions are required for hour $r\leq
t\leq T$. In our algorithm, predictions are provided by simply averaging the
last 2 weeks of data in the corresponding hour.

  * •

The hyperparameters to be learned by Algorithm 1 are the virtual electricity
price $\\{\theta_{t}\\}_{t=1,...,24}$ for 24 hours. These values are bounded
between $[0,10]$.

Optimization variables. The variables for the optimization at hour $r$
include:

  1. 1.

Net electricity grid import: $E_{t}^{\text{grid}}$, $T\geq t\geq r$

  2. 2.

Heat pump electricity usage: $E_{t}^{\text{hpC}}$, $T\geq t\geq r$

  3. 3.

Electric heater electricity usage: $E_{t}^{\text{ehH}}$, $T\geq t\geq r$

  4. 4.

Electric battery state of charge: $\text{SOC}_{t}^{\text{bat}}$, $T\geq t\geq
r$

  5. 5.

Heat storage state of charge: $\text{SOC}_{t}^{\text{H}}$, $T\geq t\geq r$

  6. 6.

Cooling storage state of charge: $\text{SOC}_{t}^{\text{C}}$, $T\geq t\geq r$

  7. 7.

Electrical storage action: $a_{t}^{\text{bat}}$, $T\geq t\geq r$

  8. 8.

Heat storage action: $a_{t}^{\text{Hsto}}$, $T\geq t\geq r$

  9. 9.

Cooling storage action: $a_{t}^{\text{Csto}}$, $T\geq t\geq r$

The actions of the policy at hour $r$ are $a_{r}^{\text{bat}}$,
$a_{r}^{\text{Hsto}}$, and $a_{r}^{\text{Csto}}$. The remaining variables are
considered auxiliary variables for planning.

Objective function. The objective function is given by:

$|E^{\text{grid}}_{t}-E^{\text{grid}}_{t-1}|+\theta_{t}E^{\text{grid}}_{t}+\sum_{t^{\prime}=t+1}^{T}\Big{(}|E^{\text{grid}}_{t^{\prime}}-E^{\text{grid}}_{t^{\prime}-1}|+\theta_{t^{\prime}}E^{\text{grid}}_{t^{\prime}}\Big{)}.$
(22)

Note that we use $e_{t}$ for $E^{\text{grid}}_{t}$ in the main text. Also, the
above objective is used in a standard minimization problem; to make it
consistent with the maximization problem in (2), we can take the negation of
the value.

Constraints. The constraints include both energy balance constraints and
technology constraints.

_Energy balance constraints:_

  * •

Electricity balance for each hour $t\geq r$:
$E_{t}^{\text{PV}}+E_{t}^{\text{grid}}=E_{t}^{\text{NS}}+E_{t}^{\text{hpC}}+E_{t}^{\text{ehH}}+a_{t}^{\text{bat}}C_{p}^{\text{bat}}$

  * •

Heat balance for each hour $t\geq r$:
$E_{t}^{\text{ehH}}=a_{t}^{\text{Hsto}}C_{p}^{\text{Hsto}}+H_{t}^{\text{bd}}$

  * •

Cooling balance for each hour $t\geq r$:
$E_{t}^{\text{hpC}}\text{COP}_{t}^{\text{C}}=a_{t}^{\text{Csto}}C_{p}^{\text{Csto}}+C_{t}^{\text{bd}}$

_Heat pump technology constraints:_

  * •

Maximum cooling for each hour $t\geq r$:
$E_{t}^{\text{hpC}}\leq E_{\text{max}}^{\text{hpC}}$

  * •

Minimum cooling for each hour $t\geq r$:
$E_{t}^{\text{hpC}}\geq 0$

_Electric heater technology constraints:_

  * •

Maximum limit for each hour $t\geq r$:
$E_{t}^{\text{ehH}}\leq E_{\text{max}}^{\text{ehH}}$

  * •

Minimum limit for each hour $t\geq r$:
$E_{t}^{\text{ehH}}\geq 0$

_Electric battery technology constraints:_

  * •

Initial SOC:
$SOC_{r}^{\text{bat}}=(1-C_{f}^{\text{bat}}SOC_{r-1}^{\text{bat}})+a_{r}^{\text{bat}}\eta^{\text{bat}}$

  * •

SOC updates for each hour $t\geq r$:
$SOC_{t}^{\text{bat}}=(1-C_{f}^{\text{bat}})SOC_{t-1}^{\text{bat}}+a_{t}^{\text{bat}}\eta^{\text{bat}}$

  * •

Action limits for each hour $t\geq r$:
$-1\leq a_{t}^{\text{bat}}\leq 1$

  * •

Bounds of SOC or each hour $t\geq r$:
$0\leq SOC_{t}^{\text{bat}}\leq 1$

_Heat storage technology constraints:_

  * •

Initial SOC:
$SOC_{r}^{\text{H}}=(1-C_{f}^{\text{Hsto}}SOC_{r-1}^{\text{H}})+a_{r}^{\text{Hsto}}\eta^{\text{Hsto}}$

  * •

SOC updates for each hour $t\geq r$:
$SOC_{t}^{\text{H}}=(1-C_{f}^{\text{Hsto}})SOC_{t-1}^{\text{H}}+a_{t}^{\text{Hsto}}\eta^{\text{Hsto}}$

  * •

Action limits or each hour $t\geq r$:
$-1\leq a_{t}^{\text{Hsto}}\leq 1$

  * •

Bounds of SOC or each hour $t\geq r$:
$0\leq SOC_{t}^{\text{H}}\leq 1$

_Cooling storage technology constraints:_

  * •

Initial SOC:
$SOC_{r}^{\text{C}}=(1-C_{f}^{\text{Csto}}SOC_{r-1}^{\text{C}})+a_{r}^{\text{Csto}}\eta^{\text{Csto}}$

  * •

SOC updates for each hour $t\geq r$:
$SOC_{t}^{\text{C}}=(1-C_{f}^{\text{Csto}})SOC_{t-1}^{\text{C}}+a_{t}^{\text{Csto}}\eta^{\text{Csto}}$

  * •

Action limits or each hour $t\geq r$:
$-1\leq a_{t}^{\text{Csto}}\leq 1$

  * •

Bounds of SOC or each hour $t\geq r$:
$0\leq SOC_{t}^{\text{C}}\leq 1$

The above optimization can be formulated as a linear program and solved
efficiently. For more implementation details, please refer to our code
(submitted as supplementary materials).

##  Appendix C Additional experimental results

###  C.1 Official results for the 2021 CityLearn Challenge

|  ZO-iRL(ours) | ICD-CA | IDLab | Breakfast Club
---|---|---|---|---
Total score | 0.944 | 1.070 | 1.070 |  1.130
total last year | 0.942 | 1.052 | 1.077 |  1.067
coord. score | 0.915 | 1.107 | 1.094 |  1.195
coord. score last year | 0.918 | 1.074 | 1.098 |  1.095
carbon emissions | 1.003 | 1.000 | 1.028 |  1.003
Table 3: Official results for the 2021 CityLearn Challenge (Nagy et al. 2021).
Here, the total score is the average of all 6 cost metrics considered in the
competition. The coordination score is the average of the first 4 metrics (see
the main paper for these metrics). Last year scores are calculated based on
the performance of the last year within the total 4-year simulation period.

###  C.2 Hyperparameters of ZO-iRL and baselines

Parameter | Value
---|---
# of parameter candidates $N_{k}$ | 3
Initial variance $\iota_{1}$ | 0.4
Guidance signal $\rho$ | specified in the main text
Guidance rate $\alpha_{k}$ | 1
Duration of one episode (hours) | 24
Range of virtual electricity price | $[0,5]$
State-action trajectory buffer size (days) | 7
Table 4: Hyperparameters for ZO-iRL. Parameter | Value
---|---
DDPG | DQN | PPO | TD3 | A2C | SAC
Learning rate | 1e-3 | 1e-4 | 3e-4 | 1e-3 | 7e-4 | 3e-4
# of epochs | NA | NA | 10 | NA | 5 | NA
Buffer size | 1e6 | 1e6 | NA | 1e6 | NA | 1e6
Batch size | 100 | 32 | 64 | 100 | NA | 256
Discount factor | 0.99 | 0.99 | 0.99 | 0.99 | 0.99 | 0.99
Table 5: Parameter values used for RL baselines. NA means not applicable. ADAM
optimizer is used for each baseline, where the policy is given by the NN
architecture with a $\tanh$ activation function and two layers of 256 units
each.

###  C.3 Results for all climate zones

Here, we provide results for all climate zones. Note that ZO-iRL performs some
random parameter exploration in the first few weeks, which results in worse
performance. However, over time, performance improves due to the guided ES, as
shown in Fig. 3.

Method | SAC | A2C | DDPG | PPO | TD3 | ZO-iRL
---|---|---|---|---|---|---
ramping cost | 1.244 (0.196) | 2.714 (0.186) | 1.327 (0.181) | 2.718 (0.079) | 1.500 (0.228) | 0.750 (0.010)
1-load factor | 1.162 (0.026) | 1.265 (0.012) | 1.192 (0.034) | 1.273 (0.022) | 1.279 (0.120) | 1.028 (0.001)
avg. daily peak | 1.195 (0.047) | 1.413 (0.024) | 1.218 (0.033) | 1.371 (0.018) | 1.351 (0.121) | 0.994 (0.003)
peak demand | 1.081 (0.001) | 1.149 (0.045) | 1.119 (0.041) | 1.419 (0.082) | 1.195 (0.161) | 0.950 (0.020)
net electric. peak | 0.992 (0.001) | 1.011 (0.002) | 0.993 (0.004) | 1.013 (0.001) | 1.099 (0.2066) | 1.008 (2e-4)
carbon emissions | 0.998 (0.001) | 1.017 (0.002) | 0.999 (0.004) | 1.023 (0.002) | 1.104 (0.203) | 1.010 (1e-4)
total score | 1.112 (0.044) | 1.428 (0.042) | 1.141 (0.047) | 1.469(0.026) | 1.255 (0.141) | 0.957 (0.004)

Table 6: Comparison of ZO-iRL and baselines for Climate Zone 2. The reported
values are the average and standard deviation (in brackets) across 10
independent runs.

Method | SAC | A2C | DDPG | PPO | TD3 | ZO-iRL
---|---|---|---|---|---|---
ramping cost | 1.098 (0.011) | 2.986 (0.223) | 1.367 (0.182) | 3.154 (0.066) | 1.263 (0.127) | 0.775 (0.004)
1-load factor | 1.140 (0.002) | 1.272 (0.009) | 1.188 (0.040) | 1.375 (0.038) | 1.240 (0.136) | 1.043 (0.001)
avg. daily peak | 1.169 (0.002) | 1.441 (0.024) | 1.242 (0.063) | 1.439 (0.025) | 1.277 (0.177) | 1.010 (0.003)
peak demand | 1.182 (0.001) | 1.272 (0.043) | 1.199 (0.022) | 1.472 (0.084) | 1.259 (0.127) | 0.952 (0.014)
net electric. peak | 0.994 (0.002) | 1.013 (0.003) | 0.996 (0.004) | 1.010 (0.001) | 1.099 (0.206) | 1.006 (3e-4)
carbon emissions | 1.000 (0.002) | 1.019 (0.003) | 1.002 (0.004) | 1.020 (0.001) | 1.104 (0.203) | 1.007 (1e-4)
total score | 1.097 (0.002) | 1.501 (0.047) | 1.166 (0.048) | 1.578(0.025) | 1.207 (0.152) | 0.966 (0.003)

Table 7: Comparison of ZO-iRL and baselines for Climate Zone 3. The reported
values are the average and standard deviation (in brackets) across 10
independent runs.

Method | SAC | A2C | DDPG | PPO | TD3 | ZO-iRL
---|---|---|---|---|---|---
ramping cost | 1.024 (0.006) | 2.814 (0.468) | 1.566 (0.230) | 3.070 (0.108) | 1.310 (0.072) | 0.739 (0.003)
1-load factor | 1.117 (0.003) | 1.229 (0.024) | 1.185 (0.036) | 1.449 (0.019) | 1.187 (0.038) | 1.013 (0.006)
avg. daily peak | 1.126 (0.002) | 1.411 (0.072) | 1.249 (0.070) | 1.429 (0.015) | 1.262 (0.051) | 1.003 (0.001)
peak demand | 1.134 (1e-4) | 1.238 (0.069) | 1.155 (0.036) | 1.444 (0.061) | 1.205 (0.099) | 0.999 (0.026)
net electric. peak | 0.987 (0.002) | 1.010 (0.006) | 0.995 (0.004) | 1.007 (0.002) | 0.990 (0.003) | 1.007 (4e-4)
carbon emissions | 0.994 (0.002) | 1.017 (0.006) | 1.000 (0.004) | 1.015 (0.001) | 0.996 (0.004) | 1.009 (8e-4)
total score | 1.064 (0.001) | 1.453 (0.103) | 1.192 (0.060) | 1.569(0.032) | 1.158 (0.042) | 0.962 (0.003)

Table 8: Comparison of ZO-iRL and baselines for Climate Zone 4. The reported
values are the average and standard deviation (in brackets) across 10
independent runs.

Method | SAC | A2C | DDPG | PPO | TD3 | ZO-iRL
---|---|---|---|---|---|---
ramping cost | 1.245 (0.162) | 2.603 (0.158) | 1.320 (0.161) | 2.673 (0.047) | 1.385 (0.144) | 0.789 (0.005)
1-load factor | 1.241 (0.055) | 1.314 (0.011) | 1.212 (0.033) | 1.332 (0.032) | 1.253 (0.043) | 1.045 (0.011)
avg. daily peak | 1.172 (0.049) | 1.347 (0.025) | 1.197 (0.043) | 1.346 (0.021) | 1.221 (0.044) | 1.004 (0.002)
peak demand | 1.285 (0.097) | 1.356 (0.021) | 1.203 (0.035) | 1.535 (0.184) | 1.231 (0.066) | 1.014 (0.022)
net electric. peak | 0.989 (0.002) | 1.003 (0.002) | 0.990 (0.006) | 1.003 (8e-4) | 0.993 (0.003) | 1.004 (0.001)
carbon emissions | 0.995 (0.002) | 1.009 (0.002) | 0.997 (0.007) | 1.015 (7e-4) | 0.999 (0.002) | 1.005 (0.001)
total score | 1.154 (0.057) | 1.438 (0.032) | 1.153 (0.041) | 1.484(0.042) | 1.180 (0.049) | 0.977 (0.005)

Table 9: Comparison of ZO-iRL and baselines for Climate Zone 5. The reported
values are the average and standard deviation (in brackets) across 10
independent runs. Figure 3: Learning curves of ZO-iRL and baselines for
Climate Zones 2–5. Note that ZO-iRL is the only method that consistently
achieves a cost below 1 across different runs in different climate zones.

###  C.4 Visualization of parameter evolution

In this section, we visualize the evolution of parameters under ZO-iRL for
some buildings. We also juxtapose the corresponding patterns of empirical peak
counts, net electricity use, electricity demand, heating demand, and cooling
demand. The empirical count is calculated for each hour as the number of times
the corresponding hour has the top 2 net electricity usage in a week. The
higher the empirical counts, the more frequent the corresponding hour has peak
usage. We also note that electricity usage is higher than electricity demand
due to the additional energy demand for heating and cooling.

(a) virtual electricity prices $\theta$

(b) net electricity usage

(c) electricity demands

(d) empirical counts of peaks

(e) heating demands

(f) cooling demands

Figure 4: Visualization of (a) parameter learning, (b) net electricity use (c)
electric loads, (d) empirical counts of peaks, (e) heating demand, and (d)
cooling demand for Building 1. It can be observed that Building 1 continues to
increase the virtual electricity prices for hours around 16–18 in response to
consistently observed peaks in those hours. Due to storage controls, the net
electricity usage pattern is smoother (spreading throughout the day) than the
demand patterns, as observed in all buildings.

(a) virtual electricity prices $\theta$

(b) net electricity usage

(c) electricity demands

(d) empirical counts of peaks

(e) cooling demands

Figure 5: Similar patterns can be observed for Building 3. Note that for this
building, there is no heating demands.

(a) virtual electricity prices $\theta$

(b) net electricity usage

(c) electricity demands

(d) empirical counts of peaks

(e) heating demands

(f) cooling demands

Figure 6: Visualization of (a) parameter learning, (b) net electricity demand,
(c) electric loads, (d) empirical counts of peaks, (e) heating demand, and (f)
cooling demand for Building 8. It can be observed that Building 8 increases
the virtual electricity price during hours 17–23 in response to high
electricity peaks. As peak issues are mitigated, virtual electricity prices
eventually decline, as can be seen after week 30.

(a) virtual electricity prices $\theta$

(b) net electricity usage

(c) electricity demands

(d) empirical counts of peaks

(e) heating demands

(f) cooling demands

Figure 7: Similar patterns can be observed for Building 6, where virtual
electricity prices rise in response to electricity peaks. For this building,
the peaks are more dispersed throughout the day.

(a) virtual electricity prices $\theta$

(b) net electricity usage

(c) electricity demands

(d) empirical counts of peaks

(e) heating demands

(f) cooling demands

Figure 8: Similar patterns are observed for Building 8.

###  C.5 Sensitivity analysis

In this section, we perform sensitivity analysis of ZO-iRL for different
parameters/design choices.

#### Initial variance of the transition distribution.

Initial variance $\iota_{1}$ determines how much randomness we inject during
each evolutionary update. Note that in each iteration, we reduce the variance
by $1/{k^{2}}$. We can observe from Table 10 that the ZO-iRL algorithm is
robust to the initial variance of the transition probability.

| initial variance
---|---
| 0.1 | 0.3 | 0.5
total score | 0.960 (5e-4) | 0.961 (7e-4) | 0.962 (0.001)
Table 10: Performance of ZO-iRL for different initial variance values.

#### Number of parameter candidates.

Here, we examine the number of candidates sampled for each update,
$N_{k}\in\\{3,5,7\\}$. In general, the performance on ZO-iRL is robust to this
parameter. There is a trade-off between the number of candidates used in each
update and the frequency of updates; as we increase $N_{k}$, we can expect to
find a better candidate in the larger pool; however, we may also decrease the
frequency of updates as the evaluation of each candidate takes one week in an
online setting. As a result, it appears that increasing the number of
candidates sampled does not help improve performance.

| number of candidates
---|---
| 3 | 5 | 7
total score | 0.959 (0.001) | 0.964 (0.002) | 0.965 (0.002)
Table 11: Performance of ZO-iRL for different numbers of sampled candidates
per update.

#### Guidance signal.

We report the sensitivity of the ZO-iRL algorithm to guidance signal
parameters. We consider variants of the guidance signal with respect to 1)
number of hours of top electricity use in the past day; 2) incremental value.
We keep the guidance learning rate fixed at $\alpha_{k}=1$. In the main text,
we consider the top-2 hours of electricity usage to be assigned a value of
0.02 (the rest hours are adjusted accordingly so that the sum over all hours
of the guidance signal is 0). Here, we consider the following variants: _(a)_
top-1 electricity hour to be assigned values of 0.02; _(b)_ top-2 electricity
hour to be assigned values of 0.04; _(c)_ top-3 electricity hour to be
assigned values of 0.04; _(d)_ top-6 electricity hour to be assigned values of
0.02. We can see from Table 12 that the proposed algorithm is robust to these
variants.

| guidance parameters
---|---
| top-1 | top-2 | top-3 | top-6
total score | 0.962 (2e-4) | 0.963 (0.003) | 0.963(0.002) | 0.975 (0.002)
Table 12: Performance of ZO-iRL for variants of guidance signals.