# Regression Equilibrium in Electricity Markets

Vladimir Dvorkin V. Dvorkin is with the Department of Electrical Engineering
and Computer Science, University of Michigan, Ann Arbor, MI. E-mail:
<EMAIL_ADDRESS>
###### Abstract

Renewable power producers participating in electricity markets build
forecasting models independently, relying on their own data, model and feature
preferences. In this paper, we argue that in renewable-dominated markets, such
an uncoordinated approach to forecasting results in substantial opportunity
costs for stochastic producers and additional operating costs for the power
system. As a solution, we introduce Regression Equilibrium—a welfare-optimal
state of electricity markets under uncertainty, where profit-seeking
stochastic producers do not benefit by unilaterally deviating from their
equilibrium forecast models. While the regression equilibrium maximizes the
private welfare, i.e., the average profit of stochastic producers across the
day-ahead and real-time markets, it also aligns with the socially optimal,
least-cost dispatch solution for the system. We base the equilibrium analysis
on the theory of variational inequalities, providing results on the existence
and uniqueness of regression equilibrium in energy-only markets. We also
devise two methods for computing the regression equilibrium: centralized
optimization and a decentralized ADMM-based algorithm that preserves the
privacy of regression datasets.

###### Index Terms:

ADMM, electricity markets, equilibrium, forecasting, renewable energy
generation, variational inequalities

## Selected Nomenclature

###  -A Optimization Parameters

$\mathbf{c}/\mathbf{C}$

$1^{\text{st}}/2^{\text{nd}}$-order generation cost coefficients

$\mathbf{c}^{+}/\mathbf{c}^{-}$

Vectors of upward$/$downward regulation cost

$\mathbf{s}/\mathbf{S}$

$1^{\text{st}}/2^{\text{nd}}$-order load shedding cost coefficients

$\widehat{\mathbf{w}}/\mathbf{w}$

Vectors of wind power forecast$/$realization

$\mathbf{d}$

Vector of nodal electricity demands

$\mathbf{F}$

Matrix of the power transfer distribution factors

$\overline{\mathbf{f}}$

Vector of maximum transmission line capacities

$\overline{\mathbf{p}}/\underline{\mathbf{p}}$

Vectors of maximum$/$minimum generation capacity

$\boldsymbol{\Gamma}$

Diagonal matrix of private revenue vs. loss weights

$\varrho$

Step size of the ADMM algorithm

$\boldsymbol{\varphi}$

Vector of wind power generation features

$\tau$

Scalar for Lasso regularization

###  -B Optimization Variables

$\mathbf{p}$

Vector of day-ahead generator dispatch

$\mathbf{r}^{+}/\mathbf{r}^{+}$

Vectors of real-time upward$/$downward regulation

$\mathbf{g}$

Stack of generator decisions, i.e.,
$\mathbf{g}=(\mathbf{p},\mathbf{r}^{+},\mathbf{r}^{+})$

$\mathbf{G}$

Generation profile, i.e., $\mathbf{G}=(\mathbf{g}_{1},\dots,\mathbf{g}_{n})$

$\boldsymbol{\ell}$

Vector of real-time load shedding

$\mathbf{L}$

Load shedding profile, i.e.,
$\mathbf{L}=(\boldsymbol{\ell}_{1},\dots,\boldsymbol{\ell}_{n})$

$\boldsymbol{\theta}$

Vector of private regression parameters

$\boldsymbol{\Theta}$

Regression profile, i.e.,
$\boldsymbol{\Theta}=(\boldsymbol{\theta}_{1},\dots,\boldsymbol{\theta}_{b})$

$\boldsymbol{\lambda}_{1}/\boldsymbol{\lambda}_{2}$

Day-ahead$/$real-time location marginal prices

###  -C Special Notation

Boldface lowercase$/$uppercase letters denote column vectors$/$matrices. Norm
$\left\|\mathbf{x}\right\|_{\mathbf{C}}^{2}=\mathbf{x}^{\top}\mathbf{C}\mathbf{x}$.
$f(x\;\\!|\;\\!y)$ is a function of $x$ with parameter $y$. Notation
$\mathbf{x}\perp\mathbf{y}$ means “orthogonal”
($\mathbf{x}\perp\mathbf{y}\Leftrightarrow\mathbf{x}^{\top}\mathbf{y}=\mathbf{0}$).
For a set $\mathbf{x}=(x_{1},\dots,x_{i},\dots,x_{n})$, $\mathbf{x}_{-i}$
collects all elements in $\mathbf{x}$ except that at position $i$. Operator
$[x]_{+}$ is the projection of $x$ onto non-negative orthant,
$\mathbin{\raisebox{0.0pt}{$\scriptstyle\otimes$}}$ is the Kronecker product,
$\mathbf{I}_{n}$ is $n\times n$ identity matrix, and vector $\mathbf{1}$
$(\mathbf{0})$ is a vector of ones (zeros).

##  I Introduction

Wholesale electricity markets are designed to maximize social welfare while
maintaining the power supply and demand in balance at all times. However, with
the increasing integration of stochastic energy resources, markets could
struggle to maintain this balance in a welfare-maximizing manner. One reason
is the lack of coordination between the day-ahead market, which is cleared
based on forecasts of renewable power generation, and the real-time market,
which offsets any forecast deviation, typically at higher costs to the system.
To improve their coordination, several mechanisms have proposed to integrate
probabilistic information on renewable generation into day-ahead markets,
e.g., using scenario-based [1, 2, 3, 4] or chance-constrained [5, 6, 7]
programming, thereby reducing the spillover from day-ahead markets with
imperfect forecasts to real-time markets. However, their adoption means
altering the existing market-clearing model and seeking a multifaceted
consensus of stakeholders (e.g., stochastic producers) on uncertainty
parameters—two currently unresolved challenges in practice.

There is potential for more seamless market coordination in enhancing
renewable power forecasts. In real-time markets with typically asymmetric
regulation costs, over- and under-predictions of stochastic generation are
penalized differently, providing wind power producers with incentives to bias
their forecasts to target cheaper regulation [8, 9, 10]. Perusing cheaper
regulation, renewable power producers implicitly achieve the least-cost
dispatch across the day-ahead and real-time markets [11]. While [8, 9, 10, 11]
focus on tuning a point or probabilistic forecasts regardless of the
underlying prediction model, the work in [12, 13, 14, 15] develops algorithms
to train machine learning (ML) models which directly map features into
decision-focused forecasts.

As electricity markets are gradually getting dominated by renewables, more
power producers strategize their forecast models. However, unlike in many
other domains for ML applications, in power systems the agent strategies are
coupled via shared grid infrastructure (physical coupling) and electricity
market clearing (economic coupling). The individual choices of prediction
models thus lead to ripple effects on the whole system and markets. The work
in [15] indeed shows how the choice of a prediction model for a single wind
farm affects locational marginal prices across the entire grid, with average
pricing errors—relative to ground truth prices—between $0.62 and $11.15 per
MWh.

#### Contributions

This paper studies how stochastic producers can cohesively strategize their ML
models in competitive electricity markets. In the setting of interest, the
model selection by each producer affects the revenues of other producers via
competitive price formation in the day-ahead and real-time markets. Since
private forecast models may be complementary or conflicting to the revenue-
seeking objectives of other participants, we use a game-theoretic lens. We
establish the existence of an equilibrium among ML models in electricity
markets, in which no market participant benefits by unilaterally deviating
from their equilibrium model parameters. Importantly, the study reveals the
role of this equilibrium in coordinating day-ahead and real-time markets and
maximizing the expected market welfare. More specifically,

  * 1.

We argue for the existence of a Nash equilibrium among private ML models that
inform trading strategies in electricity markets, particularly those applied
to wind power forecasting. We develop an equilibrium model in which profit-
seeking wind power producers strategize their forecast models in response to
day-ahead and real-time prices, converging to a profit-maximizing state with
no incentive to unilaterally change their forecast models.

  * 2.

We prove that the equilibrium regression profile—the ensemble of private
forecast models at equilibrium—not only benefits individual producers but also
leads to a socially optimal solution that minimizes the average dispatch cost.
To get this result, we reformulate the equilibrium problem as a variational
inequality (VI) problem and apply the well-established theory of VIs [16].

  * 3.

We develop two methods to compute the regression equilibrium: an equivalent
centralized optimization problem, and a decentralized algorithm based on the
alternating direction method of multipliers (ADMM) [17], in which stochastic
producers learn equilibrium forecast models locally without explicitly sharing
their training datasets.

#### Related Work

The barriers for implementing the stochastic market-clearing models in [1, 2,
3, 4, 5, 6, 7] motivated several strategies to coordinate the day-ahead and
real-time markets within existing, deterministic market-clearing models. The
strategy in [11] computes the optimal day-ahead wind power offering by
anticipating the cost of real-time re-dispatch across a fixed set of
scenarios. This strategy has shown a cost-saving potential of $\approx 25$% on
average in the New York ISO system [18]. Similarly, tuning operating reserve
[19, 20], ramping [21] and transmission capacity [22] requirements minimizes
the spillovers from day-ahead to real-time markets. More recently, [23]
proposed a data-driven strategy via cost-optimal optimization of uncertainty
sets for computing the ramping and transmission requirements in day-ahead
markets.

Importantly, the work in [11, 18, 19, 20, 21, 22, 23] offers enhancements
based on system operator’s intervention. In this paper, we establish that the
market coordination can be improved without external forces, solely when
profit-seeking stochastic producers adhere to their equilibrium ML models.

Decision-focused learning for improving the value of renewable power
predictions has gained momentum thanks to the work in [24, 13, 12, 25, 14, 23,
15]. The mainstream approaches rely on either end-to-end learning — which
leverages feedback from the downstream dispatch optimization using modified
variants of the gradient descent [13, 14, 15, 23] — or on constructing the
system-optimal loss function, subsequently used in standard training
algorithms [12, 25]. These works, however, have taken a system operator’s
perspective, and how to achieve system-optimal learning among independent
market participants remains an open question. Partially, this question was
addressed in [26], where two agents, representing two coupled infrastructure
systems, reach consensus on the single ML model for their cost-optimal
coordination. In this work, we develop algorithms to coordinate a population
of private ML learning models towards system-optimal results.

##  II Electricity Market Preliminaries

Here, we present market settings in Sec. II-A and then decision-focused wind
power forecasting in Sec. II-B.

###  II-A Market Setting

We consider a conventional market settlement with day-ahead and real-time
stages. The day-ahead market is cleared on day $d-1$ to compute the optimal
contracts for power delivery on day $d$, usually on an hourly basis. The real-
time market then settles any supply or demand deviation from the hourly day-
ahead contracts on day $d$. Similar to [11, 12], we assume:

######  Assumption 1.

Energy-only electricity market settlement.

The energy-only settlement remunerates generation units for variable costs
associated with the day-ahead dispatch and the real-time regulation; reserve,
auxiliary service and unit commitment decisions are outside the scope of the
energy-only settlement. In this convex setting, the effect of machine learning
errors on equilibrium is more evident and isolated from impacts of non-
convexities. To ease the presentation, we allow each bus to host no more than
one conventional generator, one wind power producer and one load. The loads
are modeled as fixed, price-taking demands in the day-ahead market, and as
flexible demands in the real-time market which are willing to shed their
consumption at a high cost for the system. The market welfare maximization
problem thus reduces to minimizing the social cost of electricity. We also
disregard the ramping of flexible generating units and thus drop the hour
index in the presentation. The last three assumptions are also in the spirit
of [11, 12] and are not restrictive for the main result, meaning that the
equilibrium model can be adapted to relax these assumptions in a
straightforward way. Finally, we use the matrix of power transfer distribution
factors to convert the bus power injections into power flows, thereby avoiding
explicit voltage modeling and reducing the set of equilibrium variables to
power generation. This is also non-restrictive in the DC power flow setting.

The day-ahead market-clearing then takes the form:

$\displaystyle\underset{{\mathbf{p}}}{\text{minimize}}\quad$
$\displaystyle\left\|\mathbf{p}\right\|_{\mathbf{C}}^{2}+\mathbf{c}^{\top}\mathbf{p}$
(1a) subject to
$\displaystyle\mathbf{1}^{\top}(\mathbf{p}+\widehat{\mathbf{w}}-\mathbf{d})=0,\quad\colon\mu_{1}$
(1b)
$\displaystyle|\mathbf{F}\;\\!\;\\!(\mathbf{p}+\widehat{\mathbf{w}}-\mathbf{d})|\leqslant\mathbf{\overline{f}},\quad\\!\colon\underline{\boldsymbol{\kappa}}_{1},\overline{\boldsymbol{\kappa}}_{1}$
(1c)
$\displaystyle\mathbf{\underline{p}}\leqslant\mathbf{p}\leqslant\mathbf{\overline{p}}.$
(1d)

Given the vector of electricity demand $\mathbf{d}$ and the forecast of wind
power generation $\widehat{\mathbf{w}}$, the goal here is to find the optimal
day-ahead dispatch $\mathbf{p}^{\star}$ of conventional generation, that
minimizes the quadratic generation cost $\eqref{day_ahead_mc_obj}$ satisfying
power balance (1b), power flow (1c) and generation (1d) constraints. Given the
actual realization of wind power $\mathbf{w}$, the real-time market solves the
following optimization problem:

$\displaystyle\underset{{\mathbf{r}^{+},\mathbf{r}^{-},\boldsymbol{\ell}}}{\text{minimize}}\quad$
$\displaystyle\left\|\mathbf{p}^{\star}+\mathbf{r}^{+}-\mathbf{r}^{-}\right\|_{\mathbf{C}}^{2}+\mathbf{c}^{+\top}\mathbf{r}^{+}-\mathbf{c}^{-\top}\mathbf{r}^{-}$
$\displaystyle+\left\|\boldsymbol{\ell}\right\|_{\mathbf{S}}^{2}+\mathbf{s}^{\top}\boldsymbol{\ell}$
(2a) subject to
$\displaystyle\mathbf{1}^{\top}(\mathbf{r}^{+}-\mathbf{r}^{-}+\mathbf{w}-\widehat{\mathbf{w}}+\boldsymbol{\ell})=0,\;\;\;\\!\;\;\;\;\;\\!\;\\!\colon\mu_{2}$
(2b)
$\displaystyle|\mathbf{F}(\mathbf{p}^{\star}+\widehat{\mathbf{w}}-\mathbf{d})$
$\displaystyle\\!+\mathbf{F}(\mathbf{r}^{+}-\mathbf{r}^{-}+\mathbf{w}-\widehat{\mathbf{w}}+\boldsymbol{\ell})|\leqslant\mathbf{\overline{f}},\\!\colon\\!\underline{\boldsymbol{\kappa}}_{2},\overline{\boldsymbol{\kappa}}_{2}$
(2c)
$\displaystyle\mathbf{\underline{p}}\leqslant\mathbf{p}^{\star}+\mathbf{r}^{+}-\mathbf{r}^{-}\leqslant\mathbf{\overline{p}},$
(2d)
$\displaystyle\mathbf{0}\leqslant\mathbf{r}^{+}\leqslant\mathbf{\overline{p}}-\mathbf{p}^{\star},$
(2e)
$\displaystyle\mathbf{0}\leqslant\mathbf{r}^{-}\leqslant\mathbf{p}^{\star}-\mathbf{\underline{p}},$
(2f) $\displaystyle\mathbf{0}\leqslant\boldsymbol{\ell}\leqslant\mathbf{d},$
(2g)

that minimizes the cost of re-dispatch (2a) by selecting the optimal upward
and downward regulation of generating units, $\mathbf{r}^{+}$ and
$\mathbf{r}^{-}$, respectively, and load shedding $\boldsymbol{\ell}$ in
extreme cases. An important feature of the energy-only settlement is the cost
order $\mathbf{c}^{+}\\!>\\!\mathbf{c}\\!>\\!\mathbf{c}^{-}$, which ensures
that either positive or negative deviations from the forecast always result in
additional costs to the power system. Equation (2b) requires the balance of
regulation actions w.r.t. forecast errors, and inequality (2c) guarantees that
the corresponding adjustment of day-ahead power flows remains within the line
limits. Finally, the real-time regulation is bounded by the day-ahead
contracts and technical limits via constraints (2d)–(2f), and the load
shedding is bounded by constraints (2g).

The dual variables, stated after the colon signs in problems (1) and (2), are
used for pricing day-ahead contracts and penalizing/compensating power
generation in real-time, respectively. As the two problems are strictly
monotone in their primal decision variables, their primal and dual solutions
are unique w.r.t. a particular wind power prediction $\widehat{\mathbf{w}}.$

###  II-B Revenue-Maximizing Wind Power Forecasting

Consider a single wind power producer training a machine learning model
$\mathbb{W}_{\boldsymbol{\theta}}:\mathcal{F}\mapsto\mathcal{W}$ which maps
features to power generation, i.e.,
$\mathbb{W}_{\boldsymbol{\theta}}(\boldsymbol{\varphi})=\widehat{w}$, where
$\boldsymbol{\theta}$ is a learning parameter, $\boldsymbol{\varphi}$ is the
vector of features, and $\widehat{w}$ is a prediction. The feature space
$\mathcal{F}$ includes meteorological data (wind e.g., speed and direction),
turbine data (e.g., the blade pitch angle), and their arbitrary
transformations, such as polynomial feature transformations commonly used in
the estimation of theoretical wind power curves [27]. Consider next a dataset
$\mathcal{D}=\\{(\boldsymbol{\varphi}_{1},w_{1}),\dots,(\boldsymbol{\varphi}_{n},w_{n})\\}$
of $n$ historical samples.

######  Assumption 2.

Samples size $n\gg|\boldsymbol{\varphi}|.$

Then, model $\mathbb{W}_{\boldsymbol{\theta}}$ can be optimized by minimizing
the following prediction loss function:

$\displaystyle\mathcal{L}(\boldsymbol{\theta}\;\\!|\;\\!\mathcal{D})\vcentcolon=\frac{1}{n}\sum_{i=1}^{n}\left\|\mathbb{W}_{\boldsymbol{\theta}}(\boldsymbol{\varphi}_{i})-w_{i}\right\|_{2}^{2}.$
(3)

Although the resulting model – for a given model class and data – may provide
the minimum of prediction loss, it remains myopic to regulation penalties in
the real-time market. In power systems with typically asymmetric costs, i.e.,
$|\mathbf{c}^{+}-\mathbf{c}|\neq|\mathbf{c}^{-}-\mathbf{c}|$ point-wise, wind
power producers tune their forecasts to favor less expensive regulation and
thus reduce real-time penalties [8]. The optimal forecast tuning can be
formalized as a decision-focused learning problem. Given the history of
locational marginal prices at day-ahead and real-time markets, denoted by
$\lambda_{1}$ and $\lambda_{2}$, respectively, a wind power producer optimizes
model $\mathbb{W}_{\boldsymbol{\theta}}$ by maximizing the following revenue
function:

$\displaystyle\mathcal{R}(\boldsymbol{\theta}\;\\!|\;\\!\mathcal{D})\vcentcolon=$
$\displaystyle\;\frac{1}{n}\sum_{i=1}^{n}\Big{(}\lambda_{1i}\mathbb{W}_{\boldsymbol{\theta}}(\boldsymbol{\varphi}_{i})+\lambda_{2i}(w_{i}-\mathbb{W}_{\boldsymbol{\theta}}(\boldsymbol{\varphi}_{i}))\Big{)}$
(4)

where the first summand is the day-ahead revenue from submitted forecast, and
the second summand is the real-time revenue. Thanks to the cost order
$\mathbf{c}^{+}>\mathbf{c}>\mathbf{c}^{-}$, the resulting locational marginal
prices are such that over-prediction is always penalized in real-time with
$\lambda_{2}>\lambda_{1}$, and under-prediction results in excess real-time
power supply priced at $\lambda_{2}<\lambda_{1}$. Hence, the wind power
producer has financial incentives to sell all energy at the day-ahead market,
thereby reducing real-time imbalances. However, the inherent erroneous nature
of machine learning model $\mathbb{W}_{\boldsymbol{\theta}}$ can not
completely eliminate forecast errors and thus the real-time penalties. Hence,
optimizing (4), a wind power producer finds the optimal learning parameter
$\boldsymbol{\theta}^{\star}$ that maximizes the revenue across the two
markets on average.

Both prediction- and decision-focused objectives of wind power forecasting can
be gathered in a single optimization:

$\displaystyle\underset{{\left\|\boldsymbol{\theta}\right\|_{1}\leqslant\tau}}{\text{maximize}}\quad\mathcal{R}(\boldsymbol{\theta}\;\\!|\;\\!\mathcal{D})-\gamma\mathcal{L}(\boldsymbol{\theta}\;\\!|\;\\!\mathcal{D}),$
(5)

which balances the revenue and prediction loss with some weight $\gamma$. Even
for a producer focused solely on revenue maximization, it still reasonable to
regularize the solution with a small $\gamma>0$ to prevent the model from
deviating from the fundamental physical processes underlying the data.

The variable $\boldsymbol{\theta}$ of (5) is constrained to lie in a
$\ell_{1}-$ball with $\tau>0$. This constraint forms a compact and non-empty
set for the learning model – the sufficient condition for the upcoming
equilibrium existence result. It also allows for the optimal feature selection
[28] and solution robustness [29].

##  III Regression Equilibrium

From an electricity market perspective, we are dealing with many wind power
producers optimizing their regression models on market data. We thus seek an
optimal regression profile
$\boldsymbol{\Theta}=(\boldsymbol{\theta}_{1}^{\star},\dots,\boldsymbol{\theta}_{b}^{\star})$
– a set of private learning parameters forming an equilibrium, such that no
producer would benefit from a unilateral deviation from their equilibrium
forecast model. Such regression profile defines a pure Nash equilibrium, whose
existence and computation we study in this section. We first formalize the
optimization problems of market participants and provide equilibrium
conditions in Sec. III-A. We then study equilibrium properties in Sec. III-B
and its computation in Sec. III-C.

###  III-A Regression Equilibrium Model

$\displaystyle\mathcal{P}^{\text{W}}(\boldsymbol{\Theta}\;\\!|\;\\!\boldsymbol{\lambda}_{1},\boldsymbol{\lambda}_{2},\boldsymbol{\varphi},\mathbf{w})\vcentcolon=$
$\displaystyle\;\underbrace{\boldsymbol{\lambda}_{1}^{\top}(\boldsymbol{\Theta}\boldsymbol{\varphi})}_{\text{day-
ahead
revenue}}\\!\\!\\!+\;\underbrace{\boldsymbol{\lambda}_{2}^{\top}(\mathbf{w}-\boldsymbol{\Theta}\boldsymbol{\varphi})}_{\text{real-
time
revenue}}-\underbrace{\left\|\boldsymbol{\Theta}\boldsymbol{\varphi}-\mathbf{w}\right\|_{\mathbf{\Gamma}}^{2}}_{\text{regression
loss}}$ (6a)
$\displaystyle\mathcal{P}^{\text{G}}(\mathbf{g}\;\\!|\;\\!\boldsymbol{\lambda}_{1},\boldsymbol{\lambda}_{2})\vcentcolon=$
$\displaystyle\;\underbrace{\boldsymbol{\lambda}_{1}^{\top}\mathbf{p}}_{\text{day-
ahead
revenue}}\\!\\!\\!+\;\;\;\underbrace{\boldsymbol{\lambda}_{2}^{\top}(\mathbf{r}^{+}-\mathbf{r}^{-})}_{\text{real-
time
revenue}}-\underbrace{\left(\left\|\mathbf{p}+\mathbf{r}^{+}-\mathbf{r}^{-}\right\|_{\mathbf{C}}^{2}+\mathbf{c}^{\top}\mathbf{p}+\mathbf{c}^{+\top}\mathbf{r}^{+}-\mathbf{c}^{-\top}\mathbf{r}^{-}\right)}_{\text{generation
cost}}$ (6b)
$\displaystyle\mathcal{P}^{\text{L}}(\boldsymbol{\ell}\;\\!|\;\\!\boldsymbol{\lambda}_{2})\vcentcolon=$
$\displaystyle\;\underbrace{\boldsymbol{\lambda}_{2}^{\top}\boldsymbol{\ell}}_{\text{load
shedding
payment}}\\!\\!\\!-\;\;\;\underbrace{\left(\left\|\boldsymbol{\ell}\right\|_{\mathbf{S}}+\mathbf{s}^{\top}\boldsymbol{\ell}\right)}_{\text{load
shedding cost}}$ (6c)

Throughout the paper, we make the following assumptions.

######  Assumption 3.

The class of private machine learning models
$W_{\boldsymbol{\theta}_{1}}(\boldsymbol{\varphi}),\dots,W_{\boldsymbol{\theta}_{b}}(\boldsymbol{\varphi})$
is linear in feature vector $\boldsymbol{\varphi}$.

Before proceeding, we justify this assumption. First, although the models are
linear in $\boldsymbol{\varphi}$, the non-linearity of wind power generation
can be captured by transformed features forming vector $\boldsymbol{\varphi}$.
Indeed, wind power generation is proportional to the cubed wind speed, which
can be included in the feature vector $\boldsymbol{\varphi}$, preserving the
linearity of the model. More generally, the feature vector can be composed of
non-linear kernel functions [30], e.g., Gaussian kernels. Then,
$\ell_{1}$-regularization, introduced in problem (5), can identify the set of
optimal kernels whose linear combination best explains the non-linear process
underlying the data. Second, this assumption allows for co-optimizing the
regression profile and market-clearing decisions in one optimization
problem–the property we will leverage for computing equilibrium. Other convex
models, such as input-convex neural networks [31, 32] are also of interest,
yet to be addressed in future.

######  Assumption 4.

Private models act on the same set of features.

The assumption homogenizes private models, such that market outcomes will be
explained by benefits of equilibrium and will not be obscured by heterogeneity
of models or data.

The equilibrium is formed by market participants, whose private optimization
problems are explained next, followed by the details on equilibrium
conditions.

####  III-A1 Wind Power Producers

Each producer optimizes its own machine learning model for wind power
forecasting. However, thanks to Assumptions 3 and 4 that homogenize forecast
models of all producers, we can optimize the entire regression profile
$\boldsymbol{\Theta}$ using a single optimization. The goal of this
optimization is to maximize the profit of wind power producers in response to
electricity prices from day-ahead and real-time markets. For a particular wind
power scenario, the profit function is provided in equation (7a), which is
similar to the objective function in (5) with exception that
$\mathcal{P}^{\text{W}}$ is an aggregated profit. Here, matrix
$\boldsymbol{\Gamma}=\text{daig}[\gamma_{1},\dots,\gamma_{b}]$ gathers the
private weights of all producers. Then, the profit-maximizing problem for wind
power producers is:

$\displaystyle\underset{{\boldsymbol{\Theta}=(\boldsymbol{\theta}_{1},\dots,\boldsymbol{\theta}_{b})}}{\text{maximize}}\quad$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathcal{P}^{\text{W}}(\boldsymbol{\Theta}\;\\!|\;\\!\boldsymbol{\lambda}_{1i},\boldsymbol{\lambda}_{2i},\boldsymbol{\varphi}_{i},\mathbf{w}_{i})$
(7a) subject to
$\displaystyle\boldsymbol{\Theta}\in\mathcal{O}\vcentcolon=\mathcal{O}_{1}\times\dots\times\mathcal{O}_{b},$
(7b)

that maximizes the average profit across the training datasets. The constraint
set for private regression profiles is defined as

$\displaystyle\mathcal{O}_{j}\vcentcolon=\left\\{\boldsymbol{\theta}\;\big{|}\left\|\boldsymbol{\theta}\right\|_{1}\leqslant\tau_{j},\right\\},\quad\forall
j=1,\dots,b.$ (8)

For some feature vector $\boldsymbol{\varphi}$, the profit-maximizing
forecasts are obtained from the regression profile as
$\widehat{\mathbf{w}}=\boldsymbol{\Theta}^{\star}\boldsymbol{\varphi}$.

####  III-A2 Conventional Generators

The aggregated profit function of conventional generators is provided in
equation (6b), which includes revenue streams from the day-ahead market for
nominal power supply and from the real-time market for provided regulation.
Here, vector $\mathbf{g}=(\mathbf{p},\mathbf{r}^{+},\mathbf{r}^{-})$ collects
generator decisions in both markets. For each sample $i$ in the training
dataset $\mathcal{D}$, $\mathbf{g}_{i}$ collects the corresponding dispatch
generator decisions, and we use matrix
$\mathbf{G}=(\mathbf{g}_{i},\dots,\mathbf{g}_{n})$ to collect decisions across
the training dataset. The optimization problem for conventional producers then
takes the form:

$\displaystyle\underset{{\mathbf{G}=(\mathbf{g}_{i},\dots,\mathbf{g}_{n})}}{\text{maximize}}\quad$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathcal{P}^{\text{G}}(\mathbf{g}_{i}\;\\!|\;\\!\boldsymbol{\lambda}_{1i},\boldsymbol{\lambda}_{2i})$
(9a) subject to $\displaystyle\mathbf{G}\in\mathcal{G}$ (9b)

where the constraint set $\mathcal{G}$ is defined as

$\displaystyle\mathcal{G}\vcentcolon=\left\\{\mathbf{G}\;\begin{array}[]{|l}\mathbf{\underline{p}}\leqslant\mathbf{p}_{i}+\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-}\leqslant\mathbf{\overline{p}},\\\
\mathbf{0}\leqslant\mathbf{r}_{i}^{+}\leqslant\mathbf{\overline{p}}-\mathbf{p}_{i},\\\
\mathbf{0}\leqslant\mathbf{r}_{i}^{-}\leqslant\mathbf{p}_{i}-\mathbf{\underline{p}},\\\
\forall i=1,\dots,n\end{array}\right\\}.$ (14)

####  III-A3 Loads

Whenever wind and conventional generators lack capacity to satisfy electricity
demands, the loads can shed a part of their consumption at a very high cost to
the system. The aggregated profit function for loads in this case is given in
equation (6c), and their profit-maximizing problem is

$\displaystyle\underset{{\mathbf{L}=(\boldsymbol{\ell}_{i},\dots,\boldsymbol{\ell}_{n})}}{\text{maximize}}\quad$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathcal{P}^{\text{L}}(\boldsymbol{\ell}_{i}\;\\!|\;\\!\boldsymbol{\lambda}_{2i})$
(15a) subject to $\displaystyle\mathbf{L}\in\mathcal{S}$ (15b)

where $\mathcal{S}$ is the constraint set for load shedding decisions, i.e.,

$\displaystyle\mathcal{S}\vcentcolon=\left\\{\mathbf{L}\;\begin{array}[]{|l}\mathbf{0}\leqslant\boldsymbol{\ell}_{i}\leqslant\mathbf{d},\\\
\forall i=1,\dots,n\end{array}\right\\}.$ (18)

####  III-A4 Equilibrium Conditions

These conditions build upon market settings in Sec. II-A. They couple
decision-making of market participants and yield equilibrium-supporting
prices. They first require power balance at day-ahead and real-time markets
via the following complementarity conditions:

$\displaystyle
0\leqslant\mu_{1i}\perp\mathbf{1}^{\top}(\mathbf{p}_{i}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})\geqslant
0,$ (19a) $\displaystyle
0\leqslant\mu_{2i}\perp\mathbf{1}^{\top}(\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-}-\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}+\mathbf{w}_{i}+\boldsymbol{\ell}_{i})\geqslant
0,$ (19b) where $\mu_{1i}$ and $\mu_{2i}$ are prices used to support power
balance. These conditions imply that if the power balances are satisfied with
equality, then there exist non-negative prices that incentivize sufficient
generation to maintain the system in balance. The remaining conditions are
complementarities associated with the satisfaction of network limits in the
day-ahead market
$\displaystyle\mathbf{0}\leqslant\overline{\boldsymbol{\kappa}}_{1i}\perp\mathbf{\overline{f}}-\mathbf{F}(\mathbf{p}_{i}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})\geqslant\mathbf{0},$
(19c)
$\displaystyle\mathbf{0}\leqslant\underline{\boldsymbol{\kappa}}_{1i}\perp\mathbf{\overline{f}}+\mathbf{F}(\mathbf{p}_{i}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})\geqslant\mathbf{0},$
(19d) and in the real-time market:
$\displaystyle\mathbf{0}\leqslant\overline{\boldsymbol{\kappa}}_{2i}\perp\mathbf{\overline{f}}\;\\!-\;\\!$
$\displaystyle\mathbf{F}(\mathbf{p}_{i}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})$
$\displaystyle\;\\!-\;\\!$
$\displaystyle\mathbf{F}(\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-}-\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}+\mathbf{w}_{i}+\boldsymbol{\ell}_{i})\geqslant\mathbf{0},$
(19e)
$\displaystyle\mathbf{0}\leqslant\underline{\boldsymbol{\kappa}}_{2i}\perp\mathbf{\overline{f}}\;\\!+\;\\!$
$\displaystyle\mathbf{F}(\mathbf{p}_{i}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})$
$\displaystyle\;\\!+\;\\!$
$\displaystyle\mathbf{F}(\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-}-\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}+\mathbf{w}_{i}+\boldsymbol{\ell}_{i})\geqslant\mathbf{0},$
(19f)

which similarly state that if the power flow constraints are satisfied with
equality, as in case of congestion, there exists non-negative equilibrium
prices in the corresponding market.

Dual variables
$\mu_{1i},\mu_{2i},\overline{\boldsymbol{\kappa}}_{1i},\underline{\boldsymbol{\kappa}}_{1i},\overline{\boldsymbol{\kappa}}_{2i},\underline{\boldsymbol{\kappa}}_{2i},\forall
i=1,\dots,n,$ form location marginal prices. Drawing from [33, §2], they are
derived by relating dual variables to the day-ahead and real-time energy
volumes, resulting in

$\displaystyle\boldsymbol{\lambda}_{1i}\vcentcolon=$
$\displaystyle\;\mu_{1i}\cdot\mathbf{1}-\mathbf{F}^{\top}\\!\left(\overline{\boldsymbol{\kappa}}_{1i}-\underline{\boldsymbol{\kappa}}_{1i}+\overline{\boldsymbol{\kappa}}_{2i}-\underline{\boldsymbol{\kappa}}_{2i}\right),$
(20a) $\displaystyle\boldsymbol{\lambda}_{2i}\vcentcolon=$
$\displaystyle\;\mu_{2i}\cdot\mathbf{1}-\mathbf{F}^{\top}\\!\left(\overline{\boldsymbol{\kappa}}_{2i}-\underline{\boldsymbol{\kappa}}_{2i}\right).$
(20b)

They are used in the profit functions of market participants in equations (7a)
through (6c), respectively. Importantly, as these prices come from the
complementarity conditions above, they must be seen as functions of market
participant decisions, and the regression profile
$\boldsymbol{\Theta}=(\boldsymbol{\theta}_{1},\dots,\boldsymbol{\theta}_{b})$
in particular.

####  III-A5 The Model

The Nash Equilibrium (NE) model consists of the following profit-maximizing
problems:

$\displaystyle\left\\{\begin{aligned}
&\underset{{\boldsymbol{\Theta}\in\mathcal{O}}}{\text{maximize}}\quad\frac{1}{n}\sum_{i=1}^{n}\mathcal{P}^{\text{W}}(\boldsymbol{\Theta}\;\\!|\;\\!\boldsymbol{\lambda}_{1i},\boldsymbol{\lambda}_{2i},\boldsymbol{\varphi}_{i},\mathbf{w}_{i})\\\
&\underset{{\mathbf{G}\in\mathcal{G}}}{\text{maximize}}\quad\frac{1}{n}\sum_{i=1}^{n}\mathcal{P}^{\text{G}}(\mathbf{g}_{i}\;\\!|\;\\!\boldsymbol{\lambda}_{1i},\boldsymbol{\lambda}_{2i})\\\
&\underset{{\mathbf{L}\in\mathcal{S}}}{\text{maximize}}\quad\frac{1}{n}\sum_{i=1}^{n}\mathcal{P}^{\text{L}}(\boldsymbol{\ell}_{i}\;\\!|\;\\!\boldsymbol{\lambda}_{2i})\end{aligned}\right\\}$
(21)

where the prices $\boldsymbol{\lambda}_{11},\dots,\boldsymbol{\lambda}_{1n}$
and $\boldsymbol{\lambda}_{21},\dots,\boldsymbol{\lambda}_{2n}$ are chosen
such that the following complementarities are satisfied:

$\displaystyle 0\leqslant\mu_{1i}\perp$
$\displaystyle\;\mathbf{1}^{\top}(\mathbf{p}_{i}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})\geqslant
0$ (22a) $\displaystyle 0\leqslant\mu_{2i}\perp$
$\displaystyle\;\mathbf{1}^{\top}(\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-}-\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}+\mathbf{w}_{i}+\boldsymbol{\ell}_{i})\geqslant
0$ (22b)
$\displaystyle\mathbf{0}\leqslant\overline{\boldsymbol{\kappa}}_{1i}\perp$
$\displaystyle\;\mathbf{\overline{f}}-\mathbf{F}(\mathbf{p}_{i}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})\geqslant\mathbf{0}$
(22c)
$\displaystyle\mathbf{0}\leqslant\underline{\boldsymbol{\kappa}}_{1i}\perp$
$\displaystyle\;\mathbf{\overline{f}}+\mathbf{F}(\mathbf{p}_{i}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})\geqslant\mathbf{0}$
(22d)
$\displaystyle\mathbf{0}\leqslant\overline{\boldsymbol{\kappa}}_{2i}\perp$
$\displaystyle\;\mathbf{\overline{f}}-\mathbf{F}(\mathbf{p}_{i}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})$
$\displaystyle\quad\quad\\!-\mathbf{F}(\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-}-\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}+\mathbf{w}_{i}+\boldsymbol{\ell}_{i})\geqslant\mathbf{0}$
(22e)
$\displaystyle\mathbf{0}\leqslant\underline{\boldsymbol{\kappa}}_{2i}\perp$
$\displaystyle\;\mathbf{\overline{f}}+\mathbf{F}(\mathbf{p}_{i}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})$
$\displaystyle\underbrace{\quad\quad\\!+\mathbf{F}(\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-}-\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}+\mathbf{w}_{i}+\boldsymbol{\ell}_{i})\geqslant\mathbf{0}}_{\text{shared
constraint set}\;\mathcal{E}}$ (22f)

for every sample $i=1,\dots,n$ of the training dataset $\mathcal{D}$.

###  III-B Existence and Uniqueness of Regression Equilibrium

Problem (21)-(22) is an instance of NE problem with shared constraints. To
establish the existence and uniqueness of the equilibrium, we recast it as a
variational inequality (VI) problem [16], and then rest on the well-developed
theory of VIs. The VI problem is defined as follows.

######  Definition 1.

Given a closed and convex set $\mathcal{K}\subseteq\mathbb{R}^{k}$ and a
mapping $F:\mathcal{K}\mapsto\mathbb{R}^{k}$, the VI problem, denoted
$\text{VI}(\mathcal{K},F)$, consists in finding an optimal profile
$(\boldsymbol{\Theta}^{\star},\mathbf{G}^{\star},\mathbf{L}^{\star})\in\mathcal{K}$,
called a solution to the VI problem, such that

$\left(\begin{bmatrix}[r]\boldsymbol{\Theta}\\\ \mathbf{G}\\\
\mathbf{L}\end{bmatrix}-\begin{bmatrix}[r]\boldsymbol{\Theta}^{\star}\\\
\mathbf{G}^{\star}\\\
\mathbf{L}^{\star}\end{bmatrix}\right)^{\top}F\left(\boldsymbol{\Theta}^{\star},\mathbf{G}^{\star},\mathbf{L}^{\star}\right)\geqslant
0,$

for all $(\boldsymbol{\Theta},\mathbf{G},\mathbf{L})\in\mathcal{K}$.

The NE problem (21)–(22) can recast as VI when the constraint set
$\mathcal{K}$ is an intersection of feasible regions of market participants
and the shared constraints $\mathcal{E}$.

######  Assumption 5.

Set
$\mathcal{K}\vcentcolon=\mathcal{O}\cap\mathcal{G}\cap\mathcal{S}\cap\mathcal{E}$
is compact.

This assumption simply requires that electricity demands $\mathbf{d}$, power
transfer distribution factors $\mathbf{F}$, transmission capacities
$\overline{\mathbf{f}}$, and actual wind power generation
$\mathbf{w}_{1},\dots,\mathbf{w}_{n}$ are such that there exists at least one
regression, generation and load shedding profile satisfying private and
equilibrium constraints simultaneously. The underlying mapping $F$ is defined
as

$\displaystyle
F(\boldsymbol{\Theta},\mathbf{G},\mathbf{L})=\begin{bmatrix}[l]\\\\[-8.00003pt]
\nabla_{\boldsymbol{\Theta}}\frac{1}{n}\textstyle\sum_{i=1}^{n}\mathcal{P}^{\text{W}}(\boldsymbol{\Theta}\;\\!|\;\\!\mathbf{0},\mathbf{0},\boldsymbol{\varphi}_{i},\mathbf{w}_{i})\\\\[10.00002pt]
\nabla_{\mathbf{g}_{1}}\frac{1}{n}\;\\!\mathcal{P}^{\text{G}}\\!\;\\!\;(\mathbf{g}_{1}\;\\!|\;\\!\mathbf{0},\mathbf{0})\\\\[-2.5pt]
\dots\\\\[-1.49994pt]
\nabla_{\mathbf{g}_{n}}\frac{1}{n}\mathcal{P}^{\text{G}}\\!\;\\!\;(\mathbf{g}_{n}\;\\!|\;\\!\mathbf{0},\mathbf{0})\\\\[10.00002pt]
\nabla_{\boldsymbol{\ell}_{1}}\frac{1}{n}\;\\!\mathcal{P}^{\text{L}}\\!\;\\!\;(\boldsymbol{\ell}_{1}\;\\!|\;\\!\mathbf{0},\mathbf{0})\\\\[-2.5pt]
\dots\\\\[-1.49994pt]
\nabla_{\boldsymbol{\ell}_{n}}\frac{1}{n}\mathcal{P}^{\text{L}}\\!\;\\!\;(\boldsymbol{\ell}_{n}\;\\!|\;\\!\mathbf{0},\mathbf{0})\\\
\end{bmatrix},$ (23)

i.e., it stacks the partial derivatives of the private objective functions
excluding the terms related to the shared constraints. Such a definition of
$F$ is consistent with other NE games with shared constraints, e.g., see [34,
p.46] to name one example.

Thereby, the NE problem (21)-(22) with shared constraints corresponds to the
$\text{VI}(\mathcal{K},F)$ problem defined above. This correspondence means
that if $(\boldsymbol{\Theta}^{\star},\mathbf{G}^{\star},\mathbf{L}^{\star})$
solves $\text{VI}(\mathcal{K},F),$ then there exists a set of prices
$\boldsymbol{\lambda}^{\star}=(\boldsymbol{\lambda}_{11}^{\star},\dots,\boldsymbol{\lambda}_{1n}^{\star},\boldsymbol{\lambda}_{21}^{\star},\dots,\boldsymbol{\lambda}_{2n}^{\star})$,
such that
$(\boldsymbol{\Theta}^{\star},\mathbf{G}^{\star},\mathbf{L}^{\star},\boldsymbol{\lambda}^{\star})$
is the solution to the NE problem (21)-(22). This correspondence also allows
us to establish the existence of the NE from the existence of a solution to
the VI. Furthermore, the uniqueness of the VI solution also implies the
uniqueness of the NE, as shown via the following result.

######  Theorem 1.

For the NE problem (21)–(22) resting on Assumptions 1–4, suppose that
Assumption 5 also holds. Then there exists a unique equilibrium profile
$(\boldsymbol{\Theta}^{\star},\mathbf{G}^{\star},\mathbf{S}^{\star})$, such
that for every wind power producer $j=1,\dots,b$, we have:

$\displaystyle\mathcal{P}^{\text{W}}(\boldsymbol{\theta}_{j}\;\\!|\;\\!\boldsymbol{\lambda}(\boldsymbol{\theta}_{j},\boldsymbol{\Theta}_{\\!-\\!j}^{\star},\mathbf{G}^{\star},\mathbf{S}^{\star}),\mathcal{D})$
$\displaystyle\quad\quad\quad\leqslant\mathcal{P}^{\text{W}}(\boldsymbol{\theta}_{j}^{\star}\;\\!|\;\\!\boldsymbol{\lambda}(\boldsymbol{\theta}_{j}^{\star},\boldsymbol{\Theta}_{\\!-\\!j}^{\star},\mathbf{G}^{\star},\mathbf{S}^{\star}),\mathcal{D}),\forall\boldsymbol{\theta}_{j}\in\mathcal{O}_{j}$
(24)

The proof of this result is relegated to Appendix -A. One major implication is
that there exists a unique regression equilibrium profile
$\boldsymbol{\Theta}^{\star}$ among wind power producers, such that any
unilateral deviation from the equilibrium regression would not results in
extra profits, in a competitive setting.

######  Remark 1.

The result of Theorem 1 holds on the training dataset. In Section IV, we will
demonstrate how the equilibrium regression profile generalizes on the testing
dataset as well.

###  III-C Computing Regression Equilibrium

An important upshot of Theorem 1 is the existence of a centralized
optimization problem which solves the equilibrium.

######  Corollary 1.

Under Assumptions 1 to 5, the following regularized cost-minimization problem
solves the NE (21)–(22):

$\displaystyle\underset{{\boldsymbol{\Theta},\mathbf{G},\mathbf{L}}}{\text{minimize}}\quad$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\Big{(}\left\|\mathbf{p}_{i}+\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-}\right\|_{\mathbf{C}}^{2}+\mathbf{c}^{\top}\mathbf{p}_{i}+\mathbf{c}^{+\top}\mathbf{r}_{i}^{+}$
$\displaystyle-\mathbf{c}^{-\top}\mathbf{r}_{i}^{-}+\left\|\boldsymbol{\ell}_{i}\right\|_{\mathbf{S}}^{2}+\mathbf{s}^{\top}\boldsymbol{\ell}_{i}+\left\|\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{w}_{i}\right\|_{\boldsymbol{\Gamma}}^{2}\Big{)}$
subject to
$\displaystyle\mathbf{1}^{\top}(\mathbf{p}_{i}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})=0,$
$\displaystyle\mathbf{1}^{\top}(\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-}-\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}+\mathbf{w}_{i}+\boldsymbol{\ell}_{i})=0,$
$\displaystyle|\mathbf{F}(\mathbf{p}_{i}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})|\leqslant\mathbf{\overline{f}},$
$\displaystyle|\mathbf{F}(\mathbf{p}_{i}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})$
$\displaystyle\quad+\mathbf{F}(\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-}-\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}+\mathbf{w}_{i}+\boldsymbol{\ell}_{i})|\leqslant\mathbf{\overline{f}},$
$\displaystyle\mathbf{\underline{p}}\leqslant\mathbf{p}_{i}+\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-}\leqslant\mathbf{\overline{p}},$
$\displaystyle\mathbf{0}\leqslant\mathbf{r}_{i}^{+}\leqslant\mathbf{\overline{p}}-\mathbf{p}_{i},$
$\displaystyle\mathbf{0}\leqslant\mathbf{r}_{i}^{-}\leqslant\mathbf{p}_{i}-\mathbf{\underline{p}},$
$\displaystyle\mathbf{0}\leqslant\boldsymbol{\ell}_{i}\leqslant\mathbf{d},\quad\quad\quad\;\;\;\;\\!\forall
i=1,\dots,n,$
$\displaystyle\left\|\boldsymbol{\theta}_{j}\right\|_{1}\leqslant\tau_{j},\;\quad\quad\quad\;\;\\!\forall
j=1,\dots,k.$

The objective here is to minimize the average social cost of electricity,
regularized by the total regression loss up to chosen parameter
$\boldsymbol{\Gamma}.$ The feasible region includes the day-ahead and real-
time market-clearing constraints, enforced on each sample $i$ of the training
dataset $\mathcal{D}$, as well as $\ell_{1}-$regularization constraints on
regression parameters. With a small $\boldsymbol{\Gamma}$, the optimal
regression profile $\boldsymbol{\Theta}^{\star}$ minimizes the cost across the
day-ahead and real-time markets on average, thus improving the coordination
between the day-ahead and real-time markets.

While the centralized convex optimization can be solved efficiently by
interior-point solvers, it requires all participants to gather their data—and
their training features, in particular—in a single optimization, causing
privacy risks for market participants. To preserve the privacy of training
features, we propose an ADMM-based algorithm to solve the NE problem
(21)–(22). The algorithm consists of several steps:

Step 0: Choose any $\boldsymbol{\lambda}_{1i}^{0}\geqslant 0$ and
$\boldsymbol{\lambda}_{2i}^{0}\geqslant 0$ for each sample $i=1,\dots,n$ of
the training dataset, select a step size $\varrho>0$. and set iteration index
$k\leftarrow 0$.

Step 1: If LMPs
$\boldsymbol{\lambda}_{11}^{k},\dots,\boldsymbol{\lambda}_{1n}^{k}$ and
$\boldsymbol{\lambda}_{21}^{k},\dots,\boldsymbol{\lambda}_{2n}^{k}$ satisfy a
suitable termination criteria, then stop.

Step 2: Given LMPs
$\boldsymbol{\lambda}_{11}^{k},\dots,\boldsymbol{\lambda}_{1n}^{k}$ and
$\boldsymbol{\lambda}_{21}^{k},\dots,\boldsymbol{\lambda}_{2n}^{k}$, solve the
augmented NE $(\boldsymbol{\Theta}^{k},\mathbf{G}^{k},\mathbf{L}^{k})$ with
fixed prices:

$\displaystyle\boldsymbol{\Theta}^{k}=\underset{{\boldsymbol{\Theta}\in\mathcal{O}}}{\text{argmax}}\;\frac{1}{n}\sum_{i=1}^{n}\overline{\mathcal{P}}_{\varrho}^{\text{W}}(\boldsymbol{\Theta}\;\\!|\;\\!\boldsymbol{\lambda}_{1i}^{k},\boldsymbol{\lambda}_{2i}^{k},\boldsymbol{\varphi}_{i},\mathbf{w}_{i},\mathbf{g}_{i}^{k-1}\\!,\boldsymbol{\ell}_{i}^{k-1})$
(25a)
$\displaystyle\mathbf{G}^{k}=\underset{{\mathbf{G}\in\mathcal{G}}}{\text{argmax}}\;\frac{1}{n}\sum_{i=1}^{n}\overline{\mathcal{P}}_{\varrho}^{\text{G}}(\mathbf{g}_{i}\;\\!|\;\\!\boldsymbol{\lambda}_{1i}^{k},\boldsymbol{\lambda}_{2i}^{k},\boldsymbol{\Theta}^{k-1},\boldsymbol{\ell}_{i}^{k-1})$
(25b)
$\displaystyle\mathbf{L}^{k}=\underset{{\mathbf{L}\in\mathcal{S}}}{\text{argmax}}\;\frac{1}{n}\sum_{i=1}^{n}\overline{\mathcal{P}}_{\varrho}^{\text{L}}(\boldsymbol{\ell}_{i}\;\\!|\;\\!\boldsymbol{\lambda}_{2i}^{k},,\boldsymbol{\Theta}^{k-1},\mathbf{g}_{i}^{k-1})$
(25c)

Step 3: Update the dual variables of complementarity conditions (22) for all
$i=1,\dots,n$ samples in the dataset:

$\displaystyle\mu_{1i}^{k+1}\\!\leftarrow$
$\displaystyle\\!\left[\mu_{1i}^{k}-\varrho\cdot\mathbf{1}^{\top}(\mathbf{p}_{i}^{k}+\boldsymbol{\Theta}^{k}\boldsymbol{\varphi}_{i}-\mathbf{d})\right]_{+}$
(26a) $\displaystyle\overline{\boldsymbol{\kappa}}_{1i}^{k+1}\\!\leftarrow$
$\displaystyle\\!\left[\overline{\boldsymbol{\kappa}}_{1i}^{k}-\varrho\left(\mathbf{\overline{f}}-\mathbf{F}(\mathbf{p}_{i}^{k}+\boldsymbol{\Theta}^{k}\boldsymbol{\varphi}_{i}-\mathbf{d})\right)\right]_{+}$
(26b) $\displaystyle\underline{\boldsymbol{\kappa}}_{1i}^{k+1}\\!\leftarrow$
$\displaystyle\\!\left[\underline{\boldsymbol{\kappa}}_{1i}^{k}-\varrho\left(\mathbf{\overline{f}}+\mathbf{F}(\mathbf{p}_{i}^{k}+\boldsymbol{\Theta}^{k}\boldsymbol{\varphi}_{i}-\mathbf{d})\right)\right]_{+}$
(26c) $\displaystyle\mu_{2i}^{k+1}\\!\leftarrow$
$\displaystyle\\!\left[\mu_{2i}^{k}\\!-\\!\varrho\cdot\mathbf{1}^{\top}(\mathbf{r}_{i}^{+,k}\\!\\!-\mathbf{r}_{i}^{-,k}\\!\\!-\boldsymbol{\Theta}^{k}\boldsymbol{\varphi}_{i}+\mathbf{w}_{i}+\boldsymbol{\ell}_{i}^{k})\right]_{+}$
(26d) $\displaystyle\overline{\boldsymbol{\kappa}}_{2i}^{k+1}\\!\leftarrow$
$\displaystyle\Big{[}\overline{\boldsymbol{\kappa}}_{2i}^{k}-\varrho\Big{(}\mathbf{\overline{f}}-\mathbf{F}(\mathbf{p}_{i}^{k}+\boldsymbol{\Theta}^{k}\boldsymbol{\varphi}_{i}-\mathbf{d})$
$\displaystyle\quad\quad\;\;\\!\;\\!\;\;\\!\\!\\!-\mathbf{F}(\mathbf{r}_{i}^{+,k}\\!\\!-\mathbf{r}_{i}^{-,k}\\!\\!-\boldsymbol{\Theta}^{k}\boldsymbol{\varphi}_{i}+\mathbf{w}_{i}+\boldsymbol{\ell}_{i}^{k})\Big{)}\Big{]}_{+}$
(26e) $\displaystyle\underline{\boldsymbol{\kappa}}_{2i}^{k+1}\\!\leftarrow$
$\displaystyle\Big{[}\underline{\boldsymbol{\kappa}}_{2i}^{k}-\varrho\Big{(}\mathbf{\overline{f}}+\mathbf{F}(\mathbf{p}_{i}^{k}+\boldsymbol{\Theta}^{k}\boldsymbol{\varphi}_{i}-\mathbf{d})$
$\displaystyle\quad\quad\;\;\\!\;\\!\;\;\\!\\!\\!+\mathbf{F}(\mathbf{r}_{i}^{+,k}\\!\\!-\mathbf{r}_{i}^{-,k}\\!\\!-\boldsymbol{\Theta}^{k}\boldsymbol{\varphi}_{i}+\mathbf{w}_{i}+\boldsymbol{\ell}_{i}^{k})\Big{)}\Big{]}_{+}$
(26f)

and then update electricity LMPs for all $i=1,\dots,n:$

$\displaystyle\boldsymbol{\lambda}_{1i}^{k+1}\\!\leftarrow$
$\displaystyle\mu_{1i}^{k+1}\\!\cdot\\!\mathbf{1}-\mathbf{F}^{\top}\\!\left(\overline{\boldsymbol{\kappa}}_{1i}^{k+1}\\!-\\!\underline{\boldsymbol{\kappa}}_{1i}^{k+1}\\!+\\!\overline{\boldsymbol{\kappa}}_{2i}^{k+1}\\!-\\!\underline{\boldsymbol{\kappa}}_{2i}^{k+1}\right),$
(27a) $\displaystyle\boldsymbol{\lambda}_{2i}^{k+1}\\!\leftarrow$
$\displaystyle\mu_{2i}^{k+1}\\!\cdot\\!\mathbf{1}-\mathbf{F}^{\top}\\!\left(\overline{\boldsymbol{\kappa}}_{2i}^{k+1}\\!-\\!\underline{\boldsymbol{\kappa}}_{2i}^{k+1}\right).$
(27b)

Step 4: Set $k\leftarrow k+1$ and go to Step 1.

At every iteration $k$, the algorithm computes the best response of market
participants to electricity prices in the day-ahead and real-time markets by
solving the augmented NE in (25). This problem is similar to the NE in (21)
with a difference that the private profit-maximizing optimizations are
augmented with the ADMM terms related to the feasibility of shared constraints
in (22). For example, the augmented profit function of wind power producers
takes the form:

$\displaystyle\overline{\mathcal{P}}_{\varrho}^{\text{W}}$
$\displaystyle(\boldsymbol{\Theta}\;\\!|\;\\!\boldsymbol{\lambda}_{1i},\boldsymbol{\lambda}_{2i},\boldsymbol{\varphi}_{i},\mathbf{w}_{i},\mathbf{g}_{i}^{k-1},\boldsymbol{\ell}_{i}^{k-1})\vcentcolon=$
$\displaystyle\mathcal{P}^{\text{W}}(\boldsymbol{\Theta}\;\\!|\;\\!\boldsymbol{\lambda}_{1i},\boldsymbol{\lambda}_{2i},\boldsymbol{\varphi}_{i},\mathbf{w}_{i})$
$\displaystyle\quad+\frac{\varrho}{2}\left\|\mathbf{1}^{\top}(\mathbf{p}_{i}^{k-1}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})\right\|_{2}^{2}$
$\displaystyle\quad\quad+\frac{\varrho}{2}\left\|\text{max}\left\\{\mathbf{F}(\mathbf{p}_{i}^{k-1}+\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{d})-\overline{\mathbf{f}},\mathbf{0}\right\\}\right\|_{2}^{2}$
$\displaystyle\quad\quad\quad+\dots$ (28)

where the augmented terms regularize the regression profile
$\boldsymbol{\Theta}$ to ensure the feasibility of day-ahead and real-time
market constraints. Then, the best response is used for updating the dual
variables via (26), where the current value is adjusted to account for either
power balance or flow constraint violations using a small step size $\varrho$.
These updates include a projection onto the non-negative orthant to satisfy
the complementarity conditions (22). These variables are subsequently used for
LMP updates (27), then conveyed to the next iteration. The algorithm
terminates when it meets a prescribed convergence criterion for LMPs, e.g.,
verifying whether their changes do not exceed some tolerance on average or for
each sample in the training data. The convergence of ADMM algorithms is well-
established, as discussed in [17].

Since problem (25a) is separable, it can be independently solved by each wind
power producer. Moreover, for updating the dual variables in (26), it suffices
for each wind power producer to submit only its prediction,
$\widehat{w}_{ji}=\boldsymbol{\theta}{j}^{\top}\boldsymbol{\varphi}_{i}$,
without sharing the regression parameters and feature data. Furthermore, the
system operator can carry out the price updates in equations (26)–(27), thus
eliminating the need to disclose sensitive transmission data to market
participants. Hence, the algorithm potentially preserves the privacy of the
model parameters, training features, and transmission data.

##  IV Numerical Experiments

###  IV-A Data and Settings

The wind power data is taken from [35] and includes $2\nicefrac{{1}}{{2}}$
years of readings. The readings are taken every 10 minutes and include the
active power output and corresponding features, among which the wind speed,
wind direction, and blade pitch angle are selected for forecasting. $5,000$
and $10,000$ readings are sampled for training and testing, respectively, at
random. We train a linear kernel regression with a Gaussian kernel
$\varphi(x)=\text{exp}(-\varsigma\left\|\mu-x\right\|{2}^{2})$, where
$\varsigma={0.3,0.03,0.7}$ are scale parameters for the three features
respectively, and $\mu$ is the center point. For each feature transformation,
we use ten Gaussian kernels, positioned at even intervals in the corresponding
feature domain; e.g., for wind speed, $\mu=\\{2,3,4,5,6,7,8,9,10\\}$ m$/$s. As
a result, for each wind power record, we have $|\boldsymbol{\varphi}|=30$
transformed features. To select only important feature transformations,
$\ell_{1}$ regularization with $\tau=10$ is used. We use the IEEE $24-$Bus
Reliability Test System under the active power increase (API) operating
condition [36], with a few modifications. Particularly, we install six wind
farms of $350$ MW capacity each, at buses $4,7,15,16,21,$ and $23$, covering
$38.4\%$ of the electricity demand when producing at maximum capacity. All
wind farms are granted the same dataset and features so that the private
regression models are homogeneous. Six natural gas units with a total
regulating capacity of $2.2$ GW are selected to provide regulating power,
which is sufficient to compensate for any forecast error in a feasible manner.
Generation cost are such that $\mathbf{c}^{+}=10\mathbf{c}$ and
$\mathbf{c}^{-}=0.05\mathbf{c}$, and the remaining network parameters ensure
that there is no need to shed loads, so we solely discuss the normal
operation. We compare three regression models:

  * •

Oracle: Provides perfect predictions for benchmarking forecasts and discussing
the benefits of equilibrium.

  * •

Baseline: A standard approach to forecasting, when wind power producers
minimize the prediction loss (3).

  * •

Equilibrium: Provides solution to the NE problem (21)–(22), with revenue-
seeking wind power producers optimizing their forecast models from market
data. The weights balancing the prediction loss and revenues are set uniformly
to $\gamma=10^{-4}$ for all wind power producers.

In the next section, we first present the economic benefits of equilibrium,
and then discuss computational aspects.

###  IV-B Results

Figure 1: Baseline and equilibrium forecasts as functions of wind speed at six
wind farms installed in the IEEE 24-Bus RTS.

The wind power forecasts resulting from the baseline and equilibrium
regression models are depicted in Fig. 1. Since the baseline regression solely
optimizes for prediction loss, the resulting forecasts adequately explain the
trends in data and are identical for all six wind farms. The equilibrium
forecasts, on the other hand, optimize for private welfare and may notably
deviate from the underlying data for each producer differently, as further
shown in Fig. 1. Such deviations from the ground truth data are driven by
optimization objectives and are well-known in the literature. For example, the
results of [12] and [13] show that the cost-optimal forecast consistently
under-predicts the wind power output in the day-ahead market. Here, however,
we account for network congestion, revealing that system-optimal under- or
over-prediction depends on the position of the wind farm in the power grid.
Indeed, the equilibrium regression systematically over-predicts power
generation at buses $4$ and $6$ and under-predicts power generation at buses
$7$, $17$, $21$, and $23$.

We now focus on the driving forces behind such equilibrium predictions. For
that, we first study the market revenues using two metrics: competitive ratio
(CR) and revenue incentive to deviate $\Delta\mathcal{R}$. The CR is defined
as a relative revenue performance with respect to the oracle and lies between
$0$% and $100$%. The revenue incentive to deviate $\Delta\mathcal{R}$ is
defined as the additional revenue gained from a unilateral deviation from the
current regression model. It is computed by re-solving the NE problem
(21)–(22) for each wind power producer, having the regression models of all
other wind power producers fixed, and then taking the difference. The results,
summarized in Table I, demonstrate a significant profit increase when wind
power producers select equilibrium models. Indeed, the CR increases by more
than ten percentage points for most wind power producers compared to the
baseline. The revenue incentives to deviate also reveal that the baseline
regression is highly suboptimal for each wind power producer, who can
substantially increase their revenues by a unilateral deviation from the
baseline models. The incentives to deviate from equilibrium regressions, on
the training dataset, are marginal (due to a small weight parameter $\gamma$,
used to regularize the equilibrium), asserting the result of Theorem 1. These
incentives are a bit more than marginal on the unseen testing dataset, yet
remain small compared to those in the baseline case.

Table II summarizes the prediction and system outcomes under baseline and
equilibrium regressions and contrasts them with the oracle. While the root
mean square error (RMSE) of the equilibrium regression is notably higher, the
equilibrium leads to lower dispatch costs than the baseline regression. This
asserts Corollary 1, stating that the equilibrium achieves the least-cost
dispatch across the day-ahead and real-time markets on average. Interestingly,
the equilibrium yields larger day-ahead costs, as it tends to withhold zero-
cost wind power generation from the day-ahead market. However, it also
features real-time costs that are ten times smaller than in the baseline case.
As a result, the average cost error in equilibrium is half as small as under
the baseline regression. As a side benefit, the equilibrium regression also
leads to a lower $\text{CVaR}_{10\%}$ measure of dispatch cost, defined as the
average cost across $10\%$ of the worst-case scenarios. These statistics are
similar on training and testing datasets, highlighting the ability of the
equilibrium regression to generalize beyond the training dataset.

We next study how the equilibrium drives regression model specification.
Figure 2 shows the optimal feature selection for the baseline regression,
identical for all producers, and the equilibrium selection. Although all wind
farms are granted the same set of features, the equilibrium requires them to
select features differently, depending on the wind farm’s position in the
power grid. This result highlights the importance of coordinating feature
selection to maximize private revenues, while simultaneously achieving cost-
optimal dispatch.

Table I: Average producer revenues and demand charges on the training(testing)
dataset, $

Regression | Wind power producers |  Conventional generators |  Demands
---|---|---|---
bus $4$ | bus $7$ | bus $15$ | bus $16$ | bus $21$ | bus $23$
Oracle |  $1,383$ ($1,398$) |  $913$ ($924$) |  $1,358$ ($1,373$) |  $1,369$ ($1,383$) |  $1,351$ ($1,365$) |  $1,323$ ($1,338$) |  $50,141$ ($49,967$) |  $62,729$ ($62,589$)
Baseline |  $1,081$ ($1,106$) |  $556$ ($578$) |  $1,046$ ($1,072$) |  $1,058$ ($1,084$) |  $1,038$ ($1,064$) |  $1,008$ ($1,034$) |  $52,440$ ($52,213$) |  $63,431$ ($63,308$)
CR |  $78.2$ ($79.1$) |  $60.9$ ($62.6$) |  $77$ ($78.1$) |  $77.3$ ($78.4$) |  $76.8$ ($77.9$) |  $76.2$ ($77.3$) |  $104.6$ ($104.5$) |  $101.1$ ($101.1$)
$\Delta\mathcal{R}$ |  $177$ ($173$) |  $166$ ($166$) |  $184$ ($179$) |  $185$ ($180$) |  $184$ ($179$) |  $185$ ($180$) | ——– | ——–
Equilibrium |  $1,235$ ($1,273$) |  $711$ ($722$) |  $1,208$ ($1,238$) |  $1,225$ ($1,250$) |  $1,202$ ($1,206$) |  $1,189$ ($1,186$) |  $51,876$ ($51,720$) |  $63,980$ ($63,895$)
CR |  $89.3$ ($91.1$) |  $77.9$ ($78.1$) |  $89$ ($90.2$) |  $89.5$ ($90.4$) |  $89$ ($88.4$) |  $89.9$ ($88.6$) |  $103.5$ ($103.5$) |  $102$ ($102.1$)
$\Delta\mathcal{R}$ |  $1$ ($-9$) |  $-0$ ($3$) |  $-0$ ($-8$) |  $1$ ($-7$) |  $-1$ ($52$) |  $-0$ ($56$) | ——– | ——–

CR – competitive ratio, i.e., relative performance with respect to oracle, %

$\Delta\mathcal{R}$ – average revenue incentive to unilaterally change a
regression model, $

Table II: Forecast errors, dispatch cost, and cost errors on the training(testing) dataset Regression | RMSE, MWh | Average dispatch cost, $ | Total dispatch cost error, $
---|---|---|---
total | day-ahead | real-time | average | $\text{CVaR}_{10\%}$
Oracle | ——– |  $37,246$ $(37,145)$ |  $37,246$ $(37,145)$ | ——– | ——– | ——–
Baseline |  $88$  $(88)$ |  $39,223$ $(39,071)$ |  $37,459$ $(37,373)$ |  $1,764$ $(1,698)$ |  $1,977$ $(1,925)$ |  $8,626$ $(8,120)$
Equilibrium |  $395$ $(406)$ |  $38,326$ $(38,212)$ |  $38,154$ $(38,084)$ |  $172$ $(128)$ |  $1,080$ $(1,067)$ |  $3,555$ $(3,349)$
Figure 2: The regression weights assigned to each of the 30 kernel feature
transformations under the baseline method (red) and in equilibrium (green).
Figure 3: Convergence of the wind farm profits (solid lines) to their
equilibrium values (dashed lines, see Tab. I) throughout ADMM iterations.

We conclude by discussing computational aspects. The results in this section
are obtained by solving a centralized optimization problem from Corollary 1,
which is solved using Mosek’s interior-point solver [37], requiring less than
80s to build and converge on the standard laptop. The proposed ADMM algorithm,
with the following schedule for the step size $\varrho=100_{(k\leqslant
199)},10_{(k\geqslant 200)},5_{(k\geqslant 230)},1_{(k\geqslant 275)}$, yields
identical results, as shown in Fig. 3. While this algorithm enables
equilibrium computation without disclosing private training datasets, it is
more computationally expensive than the centralized solution. Indeed, Step 2
requires solving the augmented NE (25). Each iteration takes 129.5s on
average, and the whole algorithm runs 10.8h, revealing the trade-off between
the privacy and complexity of computing the regression equilibrium. The ADMM
algorithm, however, is amenable to further improvements such as
parallelization and adaptive step size selection [17].

##  V Conclusions

This paper revealed the benefits of coordinating machine learning models in
electricity markets, particularly in the application to wind power
forecasting. It established the concept of regression equilibrium, in which
wind power producers not only maximize their profits, but also improve the
coordination of day-ahead and real-time markets, ensuring the least-cost
dispatch solution. The results on the standard IEEE test system demonstrated a
significant boost of $14$-$28\%$ in wind power profits and system cost
reduction of 3% on average (9% in the worst case) in the equilibrium state.
Moreover, although wind power generation is commonly explained by similar
features, we found that the system-optimal features differ among wind farms,
depending on their location in the grid.

###  -A Proof of Theorem 1.

Existence. Feasible regions $\mathcal{O},$ $\mathcal{G}$ and $\mathcal{S}$ are
convex, and their intersection with a convex set of equilibrium constraints
$\mathcal{E}$ is compact per Assumption 5. Moreover, the mapping $F$ is
continuous as profit functions in (6) are all differentiable. Hence, by
Corollary 2.2.5 of the existence of the solution to the VI from [16], the
equilibrium solution exists.

Uniqueness. We first show the existence of an equivalent convex optimization
problem that solves the VI, and then show the uniqueness through the
properties of this optimization. Towards finding the equivalent optimization,
we observe that the Jacobian matrix of
$F(\boldsymbol{\Theta},\mathbf{G},\mathbf{L})$, i.e.,

$\displaystyle\\!\\!\\!\nabla F(\boldsymbol{\Theta},\mathbf{G},\mathbf{L})=$
$\displaystyle\\!\\!\\!=\frac{2}{n}\begin{bmatrix}[r]\boldsymbol{\Gamma}\\!\mathbin{\raisebox{0.0pt}{$\scriptstyle\otimes$}}\\!\sum_{i=1}^{n}\boldsymbol{\varphi}_{i}\boldsymbol{\varphi}_{i}^{\top}\\!\\!\\!&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{n}\\!\mathbin{\raisebox{0.0pt}{$\scriptstyle\otimes$}}\\!\mathbf{C}&\mathbf{I}_{n}\\!\mathbin{\raisebox{0.0pt}{$\scriptstyle\otimes$}}\\!\mathbf{C}&-\mathbf{I}_{n}\\!\mathbin{\raisebox{0.0pt}{$\scriptstyle\otimes$}}\\!\mathbf{C}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{n}\\!\mathbin{\raisebox{0.0pt}{$\scriptstyle\otimes$}}\\!\mathbf{C}&\mathbf{I}_{n}\\!\mathbin{\raisebox{0.0pt}{$\scriptstyle\otimes$}}\\!\mathbf{C}&-\mathbf{I}_{n}\\!\mathbin{\raisebox{0.0pt}{$\scriptstyle\otimes$}}\\!\mathbf{C}&\mathbf{0}\\\
\mathbf{0}&-\mathbf{I}_{n}\\!\mathbin{\raisebox{0.0pt}{$\scriptstyle\otimes$}}\\!\mathbf{C}&-\mathbf{I}_{n}\\!\mathbin{\raisebox{0.0pt}{$\scriptstyle\otimes$}}\\!\mathbf{C}&\mathbf{I}_{n}\\!\mathbin{\raisebox{0.0pt}{$\scriptstyle\otimes$}}\\!\mathbf{C}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{I}_{n}\\!\mathbin{\raisebox{0.0pt}{$\scriptstyle\otimes$}}\\!\mathbf{S}\end{bmatrix}$
(29)

is a symmetric matrix. By the Symmetry Principle Theorem 1.3.1 in [16], this
Jacobian matrix renders function $F$ a gradient map, meaning that there exists
a scalar function $f(\mathbf{x})$, such that $F(\mathbf{x})=\nabla
f(\mathbf{x})$ for all $\mathbf{x}\in\mathcal{K}$. Function $f$ can be found
by integrating the gradient map as
$f(\mathbf{x})=\int_{0}^{1}F(\mathbf{x}_{0}+t(\mathbf{x}-\mathbf{x}_{0}))^{\top}(\mathbf{x}-\mathbf{x}_{0})dt$
[16, §1]. In our case, the integration of the gradient map
$F\left(\boldsymbol{\Theta},\mathbf{G},\mathbf{L}\right)$ takes the form:

$\displaystyle f(\boldsymbol{\Theta},\mathbf{G},\mathbf{L})=$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\frac{1}{n}\begin{bmatrix}[l]\\\\[-7.5pt]
2\textstyle\sum_{i=1}^{n}\gamma_{1}((t\cdot\boldsymbol{\theta}_{1})^{\top}\boldsymbol{\varphi}_{i}-w_{1i})\cdot\boldsymbol{\varphi}_{i}\\\\[-2.5pt]
\hskip 6.99997pt\dots\\\\[-0.50003pt]
2\textstyle\sum_{i=1}^{n}\gamma_{b}((t\cdot\boldsymbol{\theta}_{b})^{\top}\boldsymbol{\varphi}_{b}-w_{bi})\cdot\boldsymbol{\varphi}_{i}\\\\[5.0pt]
\mathbf{c}+2\mathbf{C}(\mathbf{p}_{1}+\mathbf{r}_{1}^{+}-\mathbf{r}_{1}^{-})\cdot
t\\\\[-2.5pt] \hskip 6.99997pt\dots\\\\[-0.50003pt]
\mathbf{c}+2\mathbf{C}(\mathbf{p}_{n}+\mathbf{r}_{n}^{+}-\mathbf{r}_{n}^{-})\cdot
t\\\\[5.0pt]
\mathbf{c}^{+}+2\mathbf{C}(\mathbf{p}_{1}+\mathbf{r}_{1}^{+}-\mathbf{r}_{1}^{-})\cdot
t\\\\[-2.5pt] \hskip 6.99997pt\dots\\\\[-0.50003pt]
\mathbf{c}^{+}+2\mathbf{C}(\mathbf{p}_{n}+\mathbf{r}_{n}^{+}-\mathbf{r}_{n}^{-})\cdot
t\\\\[5.0pt]
-\mathbf{c}^{-}-2\mathbf{C}(\mathbf{p}_{1}+\mathbf{r}_{1}^{+}-\mathbf{r}_{1}^{-})\cdot
t\\\\[-2.5pt] \hskip 6.99997pt\dots\\\\[-0.50003pt]
-\mathbf{c}^{-}-2\mathbf{C}(\mathbf{p}_{n}+\mathbf{r}_{n}^{+}-\mathbf{r}_{n}^{-})\cdot
t\\\\[5.0pt] \mathbf{s}+2\mathbf{S}\boldsymbol{\ell}_{1}\cdot t\\\\[-2.5pt]
\hskip 6.99997pt\dots\\\\[-0.50003pt]
\mathbf{s}+2\mathbf{S}\boldsymbol{\ell}_{n}\cdot t\\\\[3.99994pt]
\end{bmatrix}^{\top}\begin{bmatrix}[l]\\\\[-9.49997pt]
\boldsymbol{\theta}_{1}\\\\[-5.0pt] \hskip 1.99997pt\vdots\\\\[-1.49994pt]
\boldsymbol{\theta}_{b}\\\ \mathbf{p}_{1}\\\\[-5.0pt] \hskip
1.99997pt\vdots\\\\[-5.0pt] \mathbf{p}_{n}\\\\[1.00006pt]
\mathbf{r}_{1}^{+}\\\\[-5.0pt] \hskip 1.99997pt\vdots\\\\[-5.0pt]
\mathbf{r}_{n}^{+}\\\\[1.00006pt] \mathbf{r}_{1}^{-}\\\\[-5.0pt] \hskip
1.99997pt\vdots\\\\[-5.0pt] \mathbf{r}_{n}^{-}\\\\[1.00006pt]
\boldsymbol{\ell}_{1}\\\\[-5.0pt] \hskip 1.99997pt\vdots\\\\[-3.50006pt]
\boldsymbol{\ell}_{n}\\\ \end{bmatrix}dt=$ $\displaystyle=$
$\displaystyle\frac{1}{n}\left(\\!2\sum_{i=1}^{n}\gamma_{1}\left(\\!(\boldsymbol{\theta}_{1}^{\top}\boldsymbol{\varphi}_{i})^{\top}(\boldsymbol{\theta}_{1}^{\top}\boldsymbol{\varphi}_{i})\int_{0}^{1}tdt-
w_{1i}\boldsymbol{\theta}_{1}^{\top}\boldsymbol{\varphi}_{i}\\!\right)\\!\\!\right)$
$\displaystyle\dots$
$\displaystyle+\frac{1}{n}\left(\\!2\sum_{i=1}^{n}\gamma_{b}\\!\left(\\!\\!(\boldsymbol{\theta}_{b}^{\top}\boldsymbol{\varphi}_{i})^{\top}(\boldsymbol{\theta}_{b}^{\top}\boldsymbol{\varphi}_{i})\int_{0}^{1}tdt-
w_{1i}\boldsymbol{\theta}_{b}^{\top}\boldsymbol{\varphi}_{i}\\!\right)\\!\\!\right)$
$\displaystyle+\frac{1}{n}\Big{(}\sum_{i=1}^{n}2(\mathbf{p}_{i}+\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-})\mathbf{C}(\mathbf{p}_{i}+\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-})\int_{0}^{1}tdt$
$\displaystyle\quad\quad\quad+\mathbf{c}^{\top}\mathbf{p}_{i}+\mathbf{c}^{+\top}\mathbf{r}_{i}^{+}-\mathbf{c}^{-\top}\mathbf{r}_{i}^{-}\Big{)}$
$\displaystyle+\frac{1}{n}\Big{(}\sum_{i=1}^{n}2\boldsymbol{\ell}_{i}^{\top}\mathbf{S}\boldsymbol{\ell}_{i}\int_{0}^{1}tdt+\mathbf{s}^{\top}\boldsymbol{\ell}_{i}\Big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{n}\left(\sum_{i=1}^{n}\gamma_{1}\left((\boldsymbol{\theta}_{1}^{\top}\boldsymbol{\varphi}_{i}-w_{1i})^{2}+w_{1i}^{2}\right)\right)\hfill{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>by completing the square}}$ $\displaystyle\dots$
$\displaystyle+\frac{1}{n}\left(\sum_{i=1}^{n}\gamma_{b}\left((\boldsymbol{\theta}_{b}^{\top}\boldsymbol{\varphi}_{i}-w_{bi})^{2}+w_{bi}^{2}\right)\right)$
$\displaystyle+\frac{1}{n}\Big{(}\sum_{i=1}^{n}(\mathbf{p}_{i}+\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-})\mathbf{C}(\mathbf{p}_{i}+\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-})$
$\displaystyle\quad\quad\quad+\mathbf{c}^{\top}\mathbf{p}_{i}+\mathbf{c}^{+\top}\mathbf{r}_{i}^{+}-\mathbf{c}^{-\top}\mathbf{r}_{i}^{-}\Big{)}$
$\displaystyle+\frac{1}{n}\Big{(}\sum_{i=1}^{n}\boldsymbol{\ell}_{i}^{\top}\mathbf{S}\boldsymbol{\ell}_{i}+\mathbf{s}^{\top}\boldsymbol{\ell}_{i}\Big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\Big{(}\left\|\boldsymbol{\Theta}\boldsymbol{\varphi}_{i}-\mathbf{w}_{i}\right\|_{\boldsymbol{\Gamma}}^{2}+\mathbf{w}_{i}^{\top}\boldsymbol{\Gamma}\mathbf{w}_{i}$
$\displaystyle\quad\quad\quad\;\;+\left\|\mathbf{p}_{i}+\mathbf{r}_{i}^{+}-\mathbf{r}_{i}^{-}\right\|_{\mathbf{C}}^{2}+\mathbf{c}^{\top}\mathbf{p}_{i}$
$\displaystyle\quad\quad\quad\quad\quad\;\;+\mathbf{c}^{+\top}\mathbf{r}_{i}^{+}-\mathbf{c}^{-\top}\mathbf{r}_{i}^{-}+\left\|\boldsymbol{\ell}_{i}\right\|_{\mathbf{S}}^{2}+\mathbf{s}^{\top}\boldsymbol{\ell}_{i}\Big{)},$
(30)

which is a convex function of NE variables. Note, w.l.o.g., we used
$\mathbf{x}_{0}=\mathbf{0}$. For non-zero $\mathbf{x}_{0}$, function $f$ would
include additional terms independent from $\boldsymbol{\Theta},\mathbf{G},$
and $\mathbf{L}$. The virtue of the symmetry principle is that our
$\text{VI}(\mathcal{K},F)$ becomes a stationary point problem of the following
optimization problem:

$\displaystyle\underset{{\boldsymbol{\Theta},\mathbf{G},\mathbf{L}}}{\text{minimize}}\quad$
$\displaystyle f(\boldsymbol{\Theta},\mathbf{G},\mathbf{L})$ (31a) subject to
$\displaystyle\boldsymbol{\Theta},\mathbf{G},\mathbf{L}\in\mathcal{K},$ (31b)

such that solving problem (31) solves the VI. Since $f$ is strictly convex
($\nabla F$ in (-A) is positive-definite), the solution to problem (31) – and
hence the solution to the VI – is unique, and inequality (1) underlying the
equilibrium holds.

We finally observe that
$\mathbf{w}_{i}^{\top}\boldsymbol{\Gamma}\mathbf{w}_{i}$ in (30) is constant
and can be disregarded in optimization, so minimizing
$f(\boldsymbol{\Theta},\mathbf{G},\mathbf{L})$ is equivalent to minimizing the
average social cost of electricity, regularized by private regression losses.
Hence, the equilibrium solution aligns with the regularized social optimum.

## References

  * [1] G. Pritchard, G. Zakeri, and A. Philpott, “A single-settlement, energy-only electric power market for unpredictable and intermittent participants,” _Oper. Res._ , vol. 58, no. 4-part-2, pp. 1210–1219, 2010.
  * [2] J. M. Morales _et al._ , “Pricing electricity in pools with wind producers,” _IEEE Trans. Power Syst._ , vol. 27, no. 3, pp. 1366–1376, 2012\.
  * [3] V. M. Zavala _et al._ , “A stochastic electricity market clearing formulation with consistent pricing properties,” _Oper. Res._ , vol. 65, no. 3, pp. 557–576, 2017.
  * [4] M. Ndrio, A. N. Madavan, and S. Bose, “Pricing conditional value at risk-sensitive economic dispatch,” in _2021 IEEE Power & Energy Society General Meeting (PESGM)_. IEEE, 2021, pp. 01–05.
  * [5] Y. Dvorkin, “A chance-constrained stochastic electricity market,” _IEEE Trans. Power Syst._ , vol. 35, no. 4, pp. 2993–3003, 2019.
  * [6] R. Mieth, J. Kim, and Y. Dvorkin, “Risk-and variance-aware electricity pricing,” _Electric Power Syst. Res._ , vol. 189, p. 106804, 2020.
  * [7] A. Ratha _et al._ , “Moving from linear to conic markets for electricity,” _Eur. J. Oper. Res._ , vol. 309, no. 2, pp. 762–783, 2023\.
  * [8] G. N. Bathurst _et al._ , “Trading wind generation in short term energy markets,” _IEEE Trans. Power Syst._ , vol. 17, no. 3, pp. 782–789, 2002\.
  * [9] P. Pinson, C. Chevallier, and G. N. Kariniotakis, “Trading wind generation from short-term probabilistic forecasts of wind power,” _IEEE Trans. Power Syst._ , vol. 22, no. 3, pp. 1148–1156, 2007.
  * [10] P. Pinson, “Distributionally robust trading strategies for renewable energy producers,” _IEEE Transactions on Energy Markets, Policy and Regulation_ , vol. 1, no. 1, pp. 37–47, 2023.
  * [11] J. M. Morales _et al._ , “Electricity market clearing with improved scheduling of stochastic production,” _Eur. J. Oper. Res._ , vol. 235, no. 3, pp. 765–774, 2014.
  * [12] Y. Zhang _et al._ , “Deriving loss function for value-oriented renewable energy forecasting,” _arXiv preprint arXiv:2310.00571_ , 2023.
  * [13] ——, “Value-oriented renewable energy forecasting for coordinated energy dispatch problems at two stages,” _arXiv preprint arXiv:2309.00803_ , 2023\.
  * [14] D. Wahdany, C. Schmitt, and J. L. Cremer, “More than accuracy: end-to-end wind power forecasting that optimises the energy system,” _Electric Power Syst. Res._ , vol. 221, p. 109384, 2023.
  * [15] V. Dvorkin and F. Fioretto, “Price-aware deep learning for electricity markets,” _arXiv preprint arXiv:2308.01436_ , 2023.
  * [16] F. Facchinei and J.-S. Pang, _Finite-dimensional variational inequalities and complementarity problems_. Springer-Verlag New York, Inc., 2003.
  * [17] S. Boyd _et al._ , “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Foundations and Trends® in Machine learning_ , vol. 3, no. 1, pp. 1–122, 2011\.
  * [18] D. Zhao _et al._ , “Uncertainty-informed renewable energy scheduling: A scalable bilevel framework,” _IEEE Transactions on Energy Markets, Policy and Regulation_ , 2023.
  * [19] J. Mays, “Quasi-stochastic electricity markets,” _INFORMS Journal on Optimization_ , vol. 3, no. 4, pp. 350–372, 2021.
  * [20] V. Dvorkin, S. Delikaraoglou, and J. M. Morales, “Setting reserve requirements to approximate the efficiency of the stochastic dispatch,” _IEEE Trans. Power Syst._ , vol. 34, no. 2, pp. 1524–1536, 2018.
  * [21] B. Wang and B. F. Hobbs, “A flexible ramping product: Can it help real-time dispatch markets approach the stochastic dispatch ideal?” _Electric Power Syst. Res._ , vol. 109, pp. 128–140, 2014.
  * [22] T. V. Jensen, J. Kazempour, and P. Pinson, “Cost-optimal atcs in zonal electricity markets,” _IEEE Trans. Power Syst._ , vol. 33, no. 4, pp. 3624–3633, 2017.
  * [23] R. Mieth and H. V. Poor, “Prescribed robustness in optimal power flow,” _arXiv preprint arXiv:2310.02957_ , 2023.
  * [24] T. Carriere and G. Kariniotakis, “An integrated approach for value-oriented energy forecasting and data-driven decision-making application to renewable energy trading,” _IEEE Trans. Smart Grid_ , vol. 10, no. 6, pp. 6933–6944, 2019.
  * [25] J. Zhang, Y. Wang, and G. Hug, “Cost-oriented load forecasting,” _Electric Power Syst. Res._ , vol. 205, p. 107723, 2022.
  * [26] V. Dvorkin, “Agent coordination via contextual regression (AgentCONCUR) for data center flexibility,” _arXiv preprint arXiv:2309.16792_ , 2023.
  * [27] Y. Wang _et al._ , “Wind power curve modeling and wind power forecasting with inconsistent data,” _IEEE Transactions on Sustainable Energy_ , vol. 10, no. 1, pp. 16–25, 2018.
  * [28] R. Tibshirani, “Regression shrinkage and selection via the lasso,” _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , vol. 58, no. 1, pp. 267–288, 1996.
  * [29] D. Bertsimas and M. S. Copenhaver, “Characterization of the equivalence of robustification and regularization in linear and matrix regression,” _Eur. J. Oper. Res._ , vol. 270, no. 3, pp. 931–942, 2018.
  * [30] J. Shawe-Taylor and N. Cristianini, _Kernel methods for pattern analysis_. Cambridge university press, 2004\.
  * [31] B. Amos _et al._ , “Input convex neural networks,” in _International Conference on Machine Learning_. PMLR, 2017, pp. 146–155.
  * [32] V. Dvorkin _et al._ , “Emission-constrained optimization of gas networks: Input-convex neural network approach,” in _2023 62nd IEEE Conference on Decision and Control (CDC)_. IEEE, 2023, pp. 1575–1579.
  * [33] E. Litvinov, “Design and operation of the locational marginal prices-based electricity markets,” _IET generation, transmission & distribution_, vol. 4, no. 2, pp. 315–323, 2010.
  * [34] G. Scutari _et al._ , “Convex optimization, game theory, and variational inequality theory,” _IEEE Signal Process. Mag._ , vol. 27, no. 3, pp. 35–49, 2010.
  * [35] “Wind power forecasting,” https://www.kaggle.com/datasets/theforcecoder/wind-power-forecasting, accessed: 2024-05-19.
  * [36] S. Babaeinejadsarookolaee _et al._ , “The power grid library for benchmarking AC optimal power flow algorithms,” _arXiv preprint arXiv:1908.02788_ , 2019.
  * [37] M. ApS, “Mosek modeling cookbook,” 2020.


# Automatic off-line design of robot swarms: exploring the transferability of
control software and design methods across different platforms ††thanks: 1MK
and DGR contributed equally to this work and should be recognized as co-first
authors. Original software was implemented by MK. Chocolate and EvoStick were
developed by GF. The experiments were performed by MK with the assistance of
DGR. The manuscript was drafted by DGR and revised by MK and MB. All authors
contributed to the development of the ideas, read the manuscript, and provided
comments. The research was directed by MB. Correspondence to
<EMAIL_ADDRESS>The project has received partial funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (DEMIURGE Project, grant agreement No
681872) and from Belgium’s Wallonia-Brussels Federation through the ARC
Advanced project GbO-Guaranteed by Optimization. The authors acknowledge
support from the Belgian Fonds de la Recherche Scientifique – FNRS and from
the Colombian Ministry of Science, Technology and Innovation – Minciencias.

Miquel Kegeleirs1  IRIDIA, Université libre de Bruxelles
Brussels, Belgium
<EMAIL_ADDRESS>David Garzón Ramos1  IRIDIA, Université libre de
Bruxelles
Brussels, Belgium
<EMAIL_ADDRESS>Lorenzo Garattoni  Gianpiero Francesca  Toyota
Motor Europe
Brussels, Belgium
<EMAIL_ADDRESS>Toyota Motor Europe
Brussels, Belgium
<EMAIL_ADDRESS>Mauro Birattari ✉ IRIDIA, Université
libre de Bruxelles
Brussels, Belgium
<EMAIL_ADDRESS>
###### Abstract

Automatic off-line design is an attractive approach to implementing robot
swarms. In this approach, a designer specifies a mission for the swarm, and an
optimization process generates suitable control software for the individual
robots through computer-based simulations. Most relevant literature has
focused on effectively transferring control software from simulation to
physical robots. For the first time, we investigate (i) whether control
software generated via automatic design is transferable across robot platforms
and (ii) whether the design methods that generate such control software are
themselves transferable. We experiment with two ground mobile platforms with
equivalent capabilities. Our measure of transferability is based on the
performance drop observed when control software and/or design methods are
ported from one platform to another. Results indicate that while the control
software generated via automatic design is transferable in some cases, better
performance can be achieved when a transferable method is directly applied to
the new platform.

###### Index Terms:

Automatic design, swarm robotics, transferability, AutoMoDe, evolutionary
robotics.

##  I Introduction

A robot swarm [1, 2] is a highly redundant group of robots that operates
autonomously without relying on centralized control or external
infrastructure. Instead, swarm individuals rely on local sensing and
communication to self-organize [3]. By acting collectively, robots can
accomplish tasks that they could not accomplish individually [4].

Designing the collective behavior of a swarm is particularly challenging. No
universally applicable methodology exists for developing the control software
of the individual robots so that a desired collective behavior emerges [5].
Typically, designers manually refine control software until the desired
collective behavior is obtained. This trial-and-error design process is
costly, time-consuming, and does not guarantee reproducible or transferable
results. Automatic off-line design [6] is an appealing alternative. In this
approach, the problem of designing control software for individual robots is
re-formulated as an optimization problem. Given mission specifications and a
platform description, an optimization algorithm searches for suitable control
software for the robots. The automatic design process aims to produce control
software that maximizes swarm performance in the mission at hand, according to
a mission-specific performance metric provided as part of the mission’s formal
specification.

_Neuroevolution_ [7] is the traditional approach to the automatic off-line
design of robot swarms. In this approach, control software takes the form of
an artificial neural network with parameters tuned via artificial evolution.
More recently, _automatic modular design_ [8] has been proposed as an
alternative to neuroevolution. In modular design, an optimization algorithm
selects and tunes predefined software modules into a specific control
architecture (e.g., finite-state machines [8] or behavior trees [9]). Both in
neuroevolution and the modular approach, the design process is conducted via
computer-based simulations. The resulting control software is then transferred
to physical robots and assessed in the target environment.

Control software produced via automatic off-line design suffers from the
_reality gap_ [10]. Unavoidable differences between the design environment
(simulation) and the deployment environment (physical robots) can cause a
performance drop in the swarm [8, 11, 12, 13, 14]. When comparing design
methods, the smaller the performance drop, the greater a method’s ability to
cross the reality gap.

Building on this idea, we study the _transferability_ of control software and
automatic design methods across different robot platforms. We contend that the
process of transferring control software across different platforms gives rise
to challenges akin to those presented by the reality gap. When discussing
transferability, we consider cases where the control software is designed for
one platform and then deployed on another, or when an automatic design method
conceived for one platform is applied to another. This is reminiscent of the
reality gap problem, which emerges when designing control software in
simulation for execution on real robots.

In this preliminary study, we first investigate the conditions under which
modular and neuroevolutionary design methods conceived for a given robot
platform can produce control software for a different one. Then, we
investigate whether control software produced via automatic design can be
directly transferred from one platform to another. The two platforms
considered are ground robots endowed with equivalent sensing and actuation
capabilities.

##  II Experiments

We consider two swarms: one comprising three e-puck [15, 16] robots, and the
other, three Mercator [17] robots. We can formally describe the sensing and
actuation capabilities of the two platforms with the same reference model RM
1.2 [18]—see Table I—which defines the inputs and outputs on which control
software operates [8]. This is key to enabling the transferability between the
two platforms. The e-puck and Mercator (Figure 1) differ in size, with the
e-puck being roughly one-third the size of the Mercator. They also differ in
linear speed and sensor range. However, their speed/size and sensor-range/size
ratios are approximately the same.

Figure 1: From left to right, e-puck and Mercator robots. TABLE I: Reference
model RM 1.2. Input and output values of e-pucks’ (EP) and Mercators’ (ME)
control software.

Input |  Value–EP |  Value–ME |  Description
---|---|---|---
$\mathit{prox}$ |  $($01\text{\,}\text{/}$;$-11\text{\,}\pi$)$ |  $($01\text{\,}\text{/}$;$-11\text{\,}\pi$)$ |  proximity vector
$\mathit{light}$ |  $($01\text{\,}\text{/}$;$-11\text{\,}\pi$)$ |  $($01\text{\,}\text{/}$;$-11\text{\,}\pi$)$ |  light vector
$\mathit{gnd}$ |  $\\{\mathit{b},\mathit{g},\mathit{w}\\}$ |  $\\{\mathit{b},\mathit{g},\mathit{w}\\}$ |  ground reading
$\mathit{n}$ |  $[0,2]$ |  $[0,2]$ |  no. neighbors
$\mathit{V}$ |  $($01\text{\,}\text{/}$;$-11\text{\,}\pi$)$ |  $($01\text{\,}\text{/}$;$-11\text{\,}\pi$)$ |  neighbors vector
Output |  Value–EP |  Value–ME |  Description
$v_{k\in\\{l,r\\}}$ |  $-1010\text{\,}\mathrm{cm}\text{/}\mathrm{s}$ |  $-3030\text{\,}\mathrm{cm}\text{/}\mathrm{s}$ |  target velocity

Period of the control cycle: $0.1\text{\,}\mathrm{s}\text{/}$.
---

We design control software for e-pucks and Mercators using two automatic
methods originally conceived for the e-puck: Chocolate and EvoStick. Chocolate
[11] is a modular design method from the AutoMoDe [19] family and EvoStick [8]
is an implementation of the neuroevolutionary approach. We select these two
methods because they have been largely used in the automatic design of
collective behaviors for e-pucks [19]—both in simulation and reality. Previous
results confirm that they can generate control software for various missions,
including aggregation, foraging, and coverage [8, 11, 12, 13, 20]. In our
experiments, we consider three missions: Aggregation, Foraging, and Grid
Exploration—see Figure 2. In Aggregation, robots must aggregate on a black
spot. In Foraging, robots must iteratively travel between small black spots
and a larger white region. In Grid Exploration, robots must explore the
environment and continuously visit every cell of a grid. We adjust the size of
the environment according to the relative size of the robots: the e-pucks
operate in an environment that is one-third the size of the Mercators’ one.
The performance of the two swarms is computed using the same performance
measures. We produce a total of 120 instances of control software using
Chocolate and EvoStick—ten for each platform, mission, and method. We assess
each instance once in simulation (expected performance) and once with physical
robots (real performance). We conduct simulations in ARGoS3 [21]—a simulator
widely used in swarm robotics research.

Mercators

e-pucks

AggregationForagingGrid
Exploration$30\text{\,}\mathrm{cm}\text{/}$$100\text{\,}\mathrm{cm}\text{/}$
Figure 2: Experimental scenarios in Aggregation, Foraging, and Grid
Exploration. The workspace of the e-pucks is about one-third the size of
Mercators’ workspace. The pictures show an example of the initial
configuration at the beginning of each mission.

##  III Results and discussion

When assessed in simulation, the control software produced by both Chocolate
and EvoStick, for both e-pucks and Mercators, demonstrates meaningful behavior
and effectively performs the missions. When assessed in reality, the control
software produced by Chocolate maintains satisfactory behavior on both
platforms. In contrast, control software produced by EvoStick does not
typically reproduce simulation results on either platform. The differences
between simulation and reality are known effects of the reality gap: previous
research has shown that design methods from the AutoMoDe family are more
robust to the reality gap than those based on neuroevolution [12, 13, 14]. The
different degree of robustness to the reality gap between Chocolate and
EvoStick had been only reported for e-pucks. Here, we show that similar
results also apply to Mercators. Demonstrative videos are available for
download in the Supplementary material111Supplementary material:
https://iridia.ulb.ac.be/supp/IridiaSupp2023-001.

To investigate the extent to which Chocolate and EvoStick are transferable
from the e-puck—the platform for which they were originally conceived—to the
Mercator, we compare the performance of the control software they produce for
the two platforms. We aggregate the performance across the three missions
considered using a Friedman test [22]—see Figure 3. In simulation, the swarm
of e-pucks performs better than the one of Mercators, both for Chocolate and
EvoStick. This result was expected, as the two design methods were conceived
for the e-puck and were applied to Mercators without any adaptation. In the
experiments with the physical robots, the swarm of e-pucks performs better
than the one of Mercators when they execute control software produced by
Chocolate. On the other hand, the swarm of Mercators performs better than the
e-pucks when they execute control software produced by EvoStick. This result
was unexpected and suggests that EvoStick, although originally conceived for
the e-pucks, is more robust to the reality gap when adopted to design control
software for Mercators. This indicates that the effects of the reality gap are
not only method-dependent but also platform-dependent. By comparing the
relative performance of Chocolate and EvoStick, we can conclude that, under
the experimental conditions considered in the study, Chocolate transfers
better from e-pucks to Mercators than EvoStick—see Figure 3.

Figure 3: Comparison of control software produced by EvoStick and Chocolate
for e-pucks (EP) and Mercators (ME). We aggregate the results obtained in the
three missions using a Friedman test. For each method and platform, we present
average ranks and $95\%$ confidence intervals. The lower, the better.

We also investigate the extent to which the control software automatically
produced by Chocolate and EvoStick can be transferred between e-pucks and
Mercators. To this end, we transfer the control software produced for e-pucks
to Mercators, and _vice versa_. We assess the control software produced by the
two methods both on the platform for which it has been produced and on the
counterpart. Also in this case, we do this in simulation and reality. We
aggregate the performance across the three missions using a Friedman test—see
Figure 4. When assessed in simulation, the control software produced by
EvoStick performs better than the one produced by Chocolate. After
transferring the control software, we observe that the one produced by
Chocolate performs better than the one produced by EvoStick—see Figure 4. The
performance drop caused by transferring the control software is significantly
larger for EvoStick than for Chocolate. We observe a rank inversion between
the two design methods. It is known that EvoStick can achieve a better
performance in simulation by _overfitting_ the design process to the simulated
model of the e-puck [12, 14]. We argue that this overfitting prevented a
proper transfer of the control software to Mercator—as it happens when
transferring control software between simulation and reality. In the
experiments with the physical robots, Chocolate and EvoStick show a
performance drop after transferring the control software between e-pucks and
Mercators. However, the control software produced by Chocolate transfers
better than the one produced by EvoStick—see Figure 4.

These preliminary results indicate that automatic design methods and the
control software they produce can be transferred from one robot platform to
another—provided that the two have equivalent sensing and actuation
capabilities. The best results are obtained by transferring a method, that is,
by applying a method originally designed for a platform to another one—as
opposed to transferring the control software it produces for the original
platform to the other one. We will also investigate whether protocols to
predict the robustness of design methods to the reality gap [14] can be used
to predict the transferability of control software across platforms.

Figure 4: Comparison of control software produced by EvoStick and Chocolate
when assessed on the platform for which it has been produced (NT, for _not
transferred_) and on the counterpart (TT, for _transferred_)—regardless
whether the design was performed for e-pucks or Mercators. We aggregate the
results obtained in the three missions using a Friedman test. We present
average ranks and $95\%$ confidence intervals. The lower, the better.

## Use of AI technologies in the writing process

During the preparation of this work, the authors used OpenAI ChatGPT to
proofread and edit for language issues and stylistic inconsistencies—see
Supplementary material. After using this tool, the authors reviewed and edited
the manuscript as needed and take full responsibility for the content.

## References

  * [1] E. Şahin, “Swarm robotics: from sources of inspiration to domains of application,” in _Swarm Robotics: SAB 2004 International Workshop_ , ser. Lecture Notes in Computer Science, E. Şahin and W. M. Spears, Eds., vol. 3342. Berlin, Germany: Springer, 2005, pp. 10–20.
  * [2] G. Beni, “From swarm intelligence to swarm robotics,” in _Swarm Robotics: SAB 2004 International Workshop_ , ser. Lecture Notes in Computer Science, E. Şahin and W. M. Spears, Eds., vol. 3342. Berlin, Germany: Springer, 2005, pp. 1–9.
  * [3] M. Dorigo, M. Birattari, and M. Brambilla, “Swarm robotics,” _Scholarpedia_ , vol. 9, no. 1, p. 1463, 2014.
  * [4] M. Dorigo, G. Theraulaz, and V. Trianni, “Swarm robotics: past, present, and future [point of view],” _Proceedings of the IEEE_ , vol. 109, no. 7, pp. 1152–1165, 2021.
  * [5] M. Brambilla, E. Ferrante, M. Birattari, and M. Dorigo, “Swarm robotics: a review from the swarm engineering perspective,” _Swarm Intelligence_ , vol. 7, no. 1, pp. 1–41, 2013.
  * [6] M. Birattari, A. Ligot, D. Bozhinoski, M. Brambilla, G. Francesca, L. Garattoni, D. Garzón Ramos, K. Hasselmann, M. Kegeleirs, J. Kuckling, F. Pagnozzi, A. Roli, M. Salman, and T. Stützle, “Automatic off-line design of robot swarms: a manifesto,” _Frontiers in Robotics and AI_ , vol. 6, p. 59, 2019.
  * [7] V. Trianni, _Evolutionary Swarm Robotics_. Berlin, Germany: Springer, 2008.
  * [8] G. Francesca, M. Brambilla, A. Brutschy, V. Trianni, and M. Birattari, “AutoMoDe: a novel approach to the automatic design of control software for robot swarms,” _Swarm Intelligence_ , vol. 8, no. 2, pp. 89–112, 2014\.
  * [9] A. Ligot, J. Kuckling, D. Bozhinoski, and M. Birattari, “Automatic modular design of robot swarms using behavior trees as a control architecture,” _PeerJ Computer Science_ , vol. 6, p. e314, 2020.
  * [10] N. Jakobi, P. Husbands, and I. Harvey, “Noise and the reality gap: the use of simulation in evolutionary robotics,” in _Advances in Artificial Life: Third European Conference on Artificial Life_ , ser. Lecture Notes in Artificial Intelligence, F. Morán, A. Moreno, J. J. Merelo, and P. Chacón, Eds., vol. 929. Berlin, Germany: Springer, 1995, pp. 704–720.
  * [11] G. Francesca, M. Brambilla, A. Brutschy, L. Garattoni, R. Miletitch, G. Podevijn, A. Reina, T. Soleymani, M. Salvaro, C. Pinciroli, F. Mascia, V. Trianni, and M. Birattari, “AutoMoDe-Chocolate: automatic design of control software for robot swarms,” _Swarm Intelligence_ , vol. 9, no. 2–3, pp. 125–152, 2015.
  * [12] A. Ligot and M. Birattari, “Simulation-only experiments to mimic the effects of the reality gap in the automatic design of robot swarms,” _Swarm Intelligence_ , vol. 14, pp. 1–24, 2020.
  * [13] K. Hasselmann, A. Ligot, J. Ruddick, and M. Birattari, “Empirical assessment and comparison of neuro-evolutionary methods for the automatic off-line design of robot swarms,” _Nature Communications_ , vol. 12, p. 4345, 2021\.
  * [14] A. Ligot, “Assessing and forecasting the performance of optimization-based design methods for robot swarms: experimental protocol & pseudo-reality predictors.” Ph.D. dissertation, Université Libre de Bruxelles, Brussels, Belgium, 2023.
  * [15] F. Mondada, M. Bonani, X. Raemy, J. Pugh, C. Cianci, A. Klaptocz, S. Magnenat, J.-C. Zufferey, D. Floreano, and A. Martinoli, “The e-puck, a robot designed for education in engineering,” in _ROBOTICA 2009: Proceedings of the 9th Conference on Autonomous Robot Systems and Competitions_ , P. Gonçalves, P. Torres, and C. Alves, Eds. Castelo Branco, Portugal: Instituto Politécnico de Castelo Branco, 2009, pp. 59–65.
  * [16] L. Garattoni, G. Francesca, A. Brutschy, C. Pinciroli, and M. Birattari, “Software infrastructure for e-puck (and TAM),” IRIDIA, Université Libre de Bruxelles, Brussels, Belgium, Tech. Rep. TR/IRIDIA/2015-004, 2015.
  * [17] M. Kegeleirs, R. Todesco, D. Garzón Ramos, G. Legarda Herranz, and M. Birattari, “Mercator: hardware and software architecture for experiments in swarm SLAM,” IRIDIA, Université libre de Bruxelles, Brussels, Belgium, Tech. Rep. TR/IRIDIA/2022-012, 2022.
  * [18] K. Hasselmann, A. Ligot, G. Francesca, D. Garzón Ramos, M. Salman, J. Kuckling, F. J. Mendiburu, and M. Birattari, “Reference models for AutoMoDe,” IRIDIA, Université Libre de Bruxelles, Brussels, Belgium, Tech. Rep. TR/IRIDIA/2018-002, 2018.
  * [19] M. Birattari, A. Ligot, and G. Francesca, “AutoMoDe: a modular approach to the automatic off-line design and fine-tuning of control software for robot swarms,” in _Automated Design of Machine Learning and Search Algorithms_ , ser. Natural Computing Series, N. Pillay and R. Qu, Eds. Cham, Switzerland: Springer, 2021, pp. 73–90.
  * [20] G. Spaey, M. Kegeleirs, D. Garzón Ramos, and M. Birattari, “Evaluation of alternative exploration schemes in the automatic modular design of robot swarms,” in _Artificial Intelligence and Machine Learning: BNAIC 2019, BENELEARN 2019_ , ser. Communications in Computer and Information Science, B. Bogaerts, G. Bontempi, P. Geurts, N. Harley, B. Lebichot, T. Lenaerts, and G. Louppe, Eds. Cham, Switzerland: Springer, 2020, vol. 1196, pp. 18–33.
  * [21] C. Pinciroli, V. Trianni, R. O’Grady, G. Pini, A. Brutschy, M. Brambilla, N. Mathews, E. Ferrante, G. A. Di Caro, F. Ducatelle, M. Birattari, L. M. Gambardella, and M. Dorigo, “ARGoS: a modular, parallel, multi-engine simulator for multi-robot systems,” _Swarm Intelligence_ , vol. 6, no. 4, pp. 271–295, 2012.
  * [22] W. J. Conover, _Practical Nonparametric Statistics_ , 3rd ed., ser. Wiley Series in Probability and Statistics. New York, NY, USA: John Wiley & Sons, 1999.