# A constrained Shannon-Fano entropy coder for image storage in synthetic DNA

Xavier Pic, Marc Antonini  I3S laboratory, Côte d’Azur University and CNRS
UMR 7271, Sophia Antipolis, France
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The exponentially increasing demand for data storage has been facing more and
more challenges during the past years. The energy costs that it represents are
also increasing, and the availability of the storage hardware is not able to
follow the storage demand’s trend. The short lifespan of conventional storage
media -10 to 20 years- forces the duplication of the hardware and worsens the
situation. The majority of this storage demand concerns ”cold” data, data very
rarely accessed but that has to be kept for long periods of time. The coding
abilities of synthetic DNA, and its long durability (several hundred years),
make it a serious candidate as an alternative storage media for ”cold” data.
In this paper, we propose a variable-length coding algorithm adapted to DNA
data storage with improved performance. The proposed algorithm is based on a
modified Shannon-Fano code that respects some biochemichal constraints imposed
by the synthesis chemistry. We have inserted this code in a JPEG compression
algorithm adapted to DNA image storage [1] and we highlighted an improvement
of the compression ratio ranging from 0.5 up to 2 bits per nucleotide compared
to the state-of-the-art solution, without affecting the reconstruction
quality.

###### Index Terms:

image, compression, DNA, JPEG, entropy-coding, Shannon-Fano.

##  I Introduction

Our society is relying more and more on the digital world, and as a
consequence, the demand for digital data storage are exponentially increasing.
The duration for which this data should be preserved is also constantly
growing, mainly for legal reasons. The current tools are not sufficient for
both this reasons, first because the storage hardware cannot follow the
storage demand anymore, and secondly, because the hardware replication needed
to ensure the durability of the data storage is confronted with the same
demand problem. For those reasons, alternatives to the current data storage
technologies have to be developed quickly. Synthetic molecules, and
particularly synthetic DNA are promising candidates as alternative data
storage media [2]. Synthetic DNA molecules have huge encoding abilities (non-
synthetic DNA is encoding living beings), are very small, and can be stored
for very long periods of time. The DNA molecules are composed of a sequence of
nucleotides - adenine (A), cytosine (C), guanine (G) and thymine (T) - that
can be used to synthesize DNA and store digital data into it [3], [4].
Information theorists started to develop tools to encode data into DNA and
some of them focused on image data storage [1, 5]. Recently, the JPEG
DNA111https://jpeg.org/jpegdna/index.html research group started an
exploration on the opportunity to develop a specific image coding solution
adapted to the technologies of DNA synthesis to compress and store images into
DNA-like sequences [6]. Some compression algorithm and coders have been
developed specifically for this paradigm of data storage [7, 8, 5].

In this work we propose a variable-length quaternary encoder adapted to DNA
data storage. This encoder is based on the Shannon-Fano algorithm revisited to
work with a quaternary alphabet composed by the four letters (nucleotides)
$\\{A,C,G,T\\}$, and constrained by some biochemical restrictions. In this
paper we considered mainly the constraint on homopolymers which consists in
avoiding codewords constructed by the repetition of the same nucleotide more
than 3 times, e.g., $AAAAA$. This algorithm showed good performance compared
to one of the best state-of-the-art coder proposed in [7]. Furthermore, we
integrated the proposed quaternary coder in the image coding standard JPEG.
The modified JPEG algorithm has shown better performance when compared to the
image codec proposed in [1] which uses the Goldman coder of [7]. We finally
established bounds to identify the margin of progress still available in the
proposed variable-length coder, to identify the improvement that could benefit
the proposed DNA-based JPEG codec.

##  II Context

###  II-A DNA data storage

A workflow has been designed to store data into those DNA molecules. The first
step in this workflow is the construction of a coder able to represent data in
a sequence of symbols of the alphabet $\\{A,C,G,T\\}$. The encoded information
has to be input into a DNA molecule, through the operation of synthesis. The
molecules then are stored in hermetic capsules. To access the data, the last
operations of the workflow are used: first the DNA molecule is sequenced,
which means that we read the sequence of nucleotides from the molecule. Then
this nucleotide sequence has to be decoded to retrieve the data. In some of
those processes, (synthesis, sequencing and storage), the DNA molecules are
either manipulated or stored. During those events, a lot of factors can alter
the integrity of the molecule and introduce errors in the nucleotide sequence.
Those errors can be substitutions (a nucleotide is changed into another one,
for example an $A$ becomes a $G$) and indels (random insertions and deletions
of nucleotides in the sequence).

Furthermore, to reduce the noise level of these different operations, the DNA
coded information stream must respect some biochemical constraints on the
combinations of bases that form a DNA fragment: homopolymers, high/low $GC$
content and repeated patterns should be avoided. Designing a coding algorithm
specific to synthetic DNA data storage appeared as a solution to accurately
tackle the issues that come with the DNA storage workflow, particularly the
noise produced in the workflow. Designing compression algorithms specific to
DNA might also optimize the compression rate at which we store the data, which
is still an important topic because of the high cost of the DNA synthesis.

###  II-B Image coding into DNA: A solution based on JPEG

The codec proposed in [1] is a modified version of the classic JPEG standard
algorithm, adapted to encode data into DNA. The classic JPEG algorithm
compresses and encodes images with the use of two main encoders. Those
encoders are used to encode the DC and AC indices obtained through the
Discrete Cosine Transform (DCT) computed on the image, block by block. The
first one is the Huffman coding that encodes the run/categories of the AC
indices, and the categories of the DC indices obtained through the DCT. In
[1], this coder has been replaced by the Goldman encoder [7] that generates,
instead of binary codes, quaternary DNA-based codes. This work focuses on
improving this specific coder to increase the performance of the general
algorithm. The second coder used in the JPEG standard is a fixed-length coder
that codes the values of the DC and AC indices obtained through the DCT. In
[1], this coder has been replaced by the fixed-length quaternary coder
proposed in [8]. The figure 1 presents the whole workflow of the JPEG image
coder adapted to DNA coding that has been proposed in our previous work [1].

In the rest of the paper, we focus our work on the design of a new efficient
constrained quaternary variable-length encoder able to compete with the
algorithm proposed in [7] and used in [1].

Figure 1: Compression workflow for image coding in a quaternary representation
based on JPEG [1, 6].

##  III Huffman Tree Coding

###  III-A Background

Entropy coding methods are variable-length lossless compression methods that
represent the different symbols of the source it encodes. The codes are
estimated with the help of a frequency table. For every symbol that needs to
be encoded, the number of appearances of this symbol in the source is
computed. Finally, the length of the codes representing every symbol varies
with its probability of appearance: if a symbol is more frequent, the
associated code will be short, if it is less frequent, the code will be long.

TABLE I: Example of a binary Huffman code Symbol | Probability | Code
---|---|---
a | 0.40 | 0
b | 0.05 | 1010
c | 0.18 | 110
d | 0.07 | 1011
e | 0.20 | 111
f | 0.10 | 100

The Huffman algorithm is the optimal symbol-by-symbol entropy coding method.
Some theoretical limits on its performance have been established with the help
of Kraft’s inequality. If the source $S$ is represented by a random variable
$X$ with a probability distribution $p$, the expected length $L$ of the
Huffman code $C$ generated by this source can be defined as:

$L(C)=\sum_{x\in\Omega}p(X=x)l(C_{x})$ (1)

$C_{x}$ being the codeword representing the symbol $x$, $l(C_{x})$ the length
of this codeword and $\Omega$ the set of all possible symbols. Given $b$ the
base in which the Huffman algorithm encodes the source, the source’s Shannon
entropy $H_{b}(X)$ can be defined as:

$H_{b}(X)=\sum_{x\in X}-p(X=x)log_{b}(p(X=x))$ (2)

Shannon’s source coding theorem [9] showed that, like for any other source
code, the expected length of the $b$-ary Huffman code $C_{b}$ representing the
source $S$ of random variable $X$ follows the law:

$H_{b}(X)\leq L(C_{b})<H_{b}(X)+1$ (3)

###  III-B The Goldman coder

Nick Goldman introduced in [7] an entropy coding method adapted to DNA. This
algorithm satisfies one of the biochemical constraints presented earlier: it
doesn’t generate codes with homopolymers. Homopolymers are the repetition of
the same nucleotide base several times in a code or in a code stream (for
example $AAAAA$ or $TTTTT$). To avoid those error-prone patterns, the Goldman
code is generated by combining a ternary Huffman code and a rotating
quaternary transcoder. The ternary Huffman coding generates codes using the
alphabet {0,1,2}. The rotating transcoder translates this $\\{0,1,2\\}$
sequence into a quaternary sequence on the alphabet $\\{A,C,G,T\\}$, without
reusing the same element twice in a row. Thanks to this rotating transcoder,
no homopolymers will be generated by the Goldman coder.

TABLE II: Goldman base table for converting ternary Huffman symbols into
nucleotides

Previous nucleotide | 0 | 1 | 2
---|---|---|---
A | T | C | G
T | A | C | G
C | A | T | G
G | A | T | C

The rotating transcoder doesn’t affect the length of the codes generated by
Huffman since it transcodes the source ternary symbol per ternary symbol. For
that reason, the length of the codes for every symbol, and the expected length
of the final code, are only conditioned by the Huffman coder. Hence, we can
establish a boundary to the expected length of the entire Goldman coder
$C_{G}$:

$H_{4}(X)\leq H_{3}(X)\leq L(C_{G})$ (4)

###  III-C Limitations for constrained quaternary codes

As previously mentioned, the theory has shown that the Huffman tree coding
algorithm was the best entropy coding algorithm when source symbols are
independent and identically distributed, i.i.d. source. This means that the
quaternary Huffman tree coding algorithm will give the best performance in
quaternary entropy coding of a i.i.d. source. The algorithm in itself is the
theoretical limit for any quaternary entropy coding, hence the quaternary
Huffman code $C_{H4}$ is a bound for any quaternary entropy coder $C$:

$\forall C\in Ce_{4},H_{4}(X)\leq L(C_{H4})\leq L(C)$ (5)

$Ce_{4}$ being the set of all the possible quaternary entropy coders. The
problem for DNA data storage is that, as we said previously, an unconstrained
code will generate a lot of noise. This was the reasoning behind the
development of the Goldman algorithm, and is still a problem that has to be
considered. Thus, one cannot simply use the quaternary Huffman coder. This
means that to develop a new entropy coder for DNA with improved performance
over the Goldman algorithm or any ternary entropy coder $C$, a good objective
in terms of performance would be to design a code $C$ such as:

$L(C_{H4})\leq L(C)\leq H_{3}(X)$ (6)

##  IV Proposed variable length quaternary encoder

###  IV-A Motivations

The goal of the Goldman coder [7] is to constrain the code and avoid the
apparition of homopolymers. The experiments have shown that the biochemical
processes of DNA data storage (synthesis, storage and sequencing) can tolerate
homopolymers until a certain length (3 for Illumina sequencers and around 5
for Nanopore sequencers), otherwise, the noise starts increasing dramatically
[10]. This means that the Goldman coder, that does not allow two consecutive
identical nucleotides, is too restrictive when avoiding homopolymers.

In order to meet the homopolymer constraint without losing too much coding
potential, we proposed a solution that consists in a quaternary tree-coder
constrained at certain positions to avoid too long homopolymers. Since tree
structures are used for the code construction, the most simple solution was to
constraint the arity (the number of sons of each node) of the tree at certain
depths. The classic Huffman algorithm builds the coding tree from the leaves
up to the root, and does not hold any information concerning the final depth
of the tree nor the final depth of any of the nodes in the tree. Designing our
algorithm from the Huffman algorithm, where there is a local control of the
tree structure, did not appear as the most efficient solution, especially for
obtaining a balanced tree.
Instead, we based our solution on the classic Shannon-Fano recursive algorithm
which builds the coding tree from the root, meaning that the current depth can
be propagated through the recursions.

###  IV-B Proposed constrained tree coder

The coder is implemented with two algorithms. The algorithm 1 is a simple
initialization: it starts the whole construction of the tree, at the root (at
depth 1). The algorithm 2 manages the recursive calls in function of two
parameters: a global one, the maximum length of homopolymer that is allowed
$max_{hl}$, and a local one, the depth $d$ at which the call was made. If the
depth $d$ is in the class 0 of congruence modulo $max_{hl}$, the tree will be
constrained to ternary, otherwise, it will not be constrained. The figure 2
shows an example of such a coding tree.

Figure 2: A section of a coding tree generated by the proposed algorithm for
$max_{h}l=3$, for readability, two of the top level branches and most branches
at depth 4 to 5 have been deleted.

In the algorithms 1 and 2, we assume that the symbols in the set $S$ and their
frequency table $f$ have both already been sorted by decreasing frequency. In
the algorithm 2, the function $partition(S,f,n)$ slices into $n$ divisions the
set of symbols $S$ and frequencies $f$, with the cumulated frequencies of
every slice being as close as possible to each other. The function
$merge\\_trees(P)$ takes the list $P$ containing 3 or 4 trees and generates a
tree with a new root, and each son of this root being one of the trees in the
input sequence $P$. The first son is labeled with a 0, the second with a 1,
the third with 2 and the last one with a 3.

return $Sfc\\_aux(S,f,max_{hl},1)$

Algorithm 1 $Sfc(S,f,max_{hl})$: Create the DNA coding tree $T$ for a set of
symbols $S$, $f$ their frequency table and $max_{hl}$ the maximum length of
homopolymer that is allowed.

if $length(S)=0$ then

return $\emptyset$

else if $length(S)=1$ then

return $create\\_leaf(S)$

else if $d\equiv 0\pmod{max_{hl}}$ then

$PS,Pf\leftarrow partition(S,f,3)$ //Constrained case

else

$PS,Pf\leftarrow partition(S,f,4)$ //Unconstrained case

end if

$PT\leftarrow$ empty list of length $length(PS)$

for $i\leftarrow 1$ to $length(PS)$ do

$PT[i]\leftarrow Sfc\\_aux(PS[i],Pf[i],max_{hl},d+1)$

end for

$T\leftarrow merge\\_trees(PT)$

return $T$

Algorithm 2 $Sfc\\_aux(S,f,max_{hl},d)$: Create the DNA coding tree at depth
$d$.

###  IV-C Transcoding into DNA

The algorithm described earlier generates a code containing for every symbol a
sequence of bases: 0, 1, 2, and 3. These sequences still need to be translated
into DNA sequences respecting the homopolymer constraint described earlier.

To transcode the sequences into DNA, two different kinds of transcoding tables
are used. According to the tree structure, the $k^{th}$ symbol in the sequence
is a ternary symbol if $k\equiv 0\pmod{max_{hl}}$ and a quaternary one
otherwise, $max_{hl}$ being the maximum homopolymer length. When the symbol is
ternary, it is translated into a nucleotide with the Goldman transcoding table
(Table II), otherwise, the quaternary trancoding table (Table III) is used.

TABLE III: Quaternary transcoding table Symbol | 0 | 1 | 2 | 3
---|---|---|---|---
Nucleotide | A | T | C | G

###  IV-D Property of the proposed tree coder

If $max_{hl}$ is the maximum homopolymer length that can be accepted,
constraining the coding tree to ternary at every depth that is a multiple of
$max_{hl}$, would prevent the apparition of homopolymers of length more than
$max_{hl}$ (see Property 1).

######  Proposition 1.

Let $C_{SFC}$ be a code obtained thanks to the proposed coder (algorithm 1).
Let $\mathcal{S}$$(C_{x})$ be the set of homopolymer patterns found in the
codeword $C_{x}$ of the symbol $x$. Let $\omega_{h}\in\mathcal{S}(C_{x})$ be a
homopolymer subsection of $(C_{x})$. Let $d_{i}$, and respectively $d_{j}$ be
the depth at wich the section $\omega_{h}$ starts, respectively ends, in the
coding tree.
Then, $\exists!{k}\in{\mathbf{N}}$ with,

$\begin{split}k\times max_{hl}&<d_{i}<(k+1)\times max_{hl}\\\ k\times
max_{hl}&<d_{j}\leq(k+1)\times max_{hl}\end{split}$ (7)

###### Proof.

By contradiction: If a code $C_{SFC}$ has a codeword $C_{x}$ containing a
homopolymer pattern $\omega_{h}$ of length bigger than $max_{hl}$, the
homopolymer can be located in the coding tree between depths $d_{i}$ and
$d_{j}$ with $l(\omega_{h})=d_{j}-d_{i}>max_{hl}$. By construction of the
coding tree, there is a constrained depth in the middle of the homopolymer
pattern, where the Goldman transcoding table II is used, hence two consecutive
nucleotides are not identical, so $\omega_{h}$ is not a homopolymer pattern. ∎

##  V Experimental results

###  V-A Experiments on synthetic sources

Before integrating the coder inside a complete image codec, we estimate the
raw gain that we can obtain in comparison to the Goldman coder ($C_{G}$). We
also compare our results to a Huffman-based constrained quaternary code
($C_{H4C}$) that we implemented, even if as mentionned in the motivations
subsection (IV-A) it was not considered as a good candidate from the
beginning.

The metric used to compare the performance of variable-length coders is the
expected length, that has been previously introduced in equation (1).

We used i.i.d Gaussian sources for those experiments. It was made sure to
generate enough samples in the sources to approach the Gaussian distribution:
100 relisations, each one containing 10000 samples. Each sample is quantized
using a dictionnary composed of 162 symbols. Since the Goldman coder is based
on a Huffman ternary tree coder, comparing the expected length of the two
coders to the ternary entropy (equation (2) with $b=3$) is appropriate: if the
new coder beats the ternary entropy, then it is better than any tenary-based
entropy coder for DNA (see equation (6)).

TABLE IV: Length expectancy results of the entropy coders, expressed in nucleotides per symbol, on Gaussian sources, with $max_{hl}=3$. $H_{4}(X)$ | $L(C_{H4})$ | $\boldsymbol{L(C_{SFC})}$ | $L(C_{H4C})$ | $H_{3}(X)$ | $L(C_{G})$
---|---|---|---|---|---
3.48 | 3.56 | 3.81 | 4.31 | 4.39 | 4.45

The length expectancy presented in table IV shows a gain of 0.64 nucleotide
per symbol for the proposed algorithm ($C_{SFC}$) over Goldman ($C_{G}$), and
a difference of 0.25 nucleotide per symbol with the unconstrained quaternary
Huffman coder ($C_{H4}$). As expected in the motivations subsection (IV-A), in
comparison, the constrained Huffman-based coder $C_{H4C}$ clearly
underperforms. Not having a control over the depth of the coding tree means
that it is harder to balance it and optimize the expected length.

As a last experiment, the performance of the coder was evaluated on a source
following a distribution that is similar to the distribution of the AC indices
produced by a JPEG codec in most natural images. The results presented in the
table V have been generated using a source that follows a frequency table of
the AC coefficients computed on the Kodak image
dataset222http://r0k.us/graphics/kodak/. Here again, the proposed constrained
Shannon-Fano quaternary encoder is the best performing.

TABLE V: Length expectancy results of the entropy coders, expressed in nucleotides per symbol, on sources using AC values frequencies, with $max_{hl}=3$ $H_{4}(X)$ | $L(C_{H4})$ | $\boldsymbol{L(C_{SFC})}$ | $H_{3}(X)$ | $L(C_{H4C})$ | $L(C_{G})$
---|---|---|---|---|---
1.21 | 1.39 | 1.43 | 1.53 | 1.61 | 1.64

###  V-B Experiments on image coding

As explained in the introduction, the variable length coder is only one of the
coders used in the JPEG DNA experimentation software’s algorithm [6]. It is
used to encode the category of a DC coefficient or the run/category (a
combination of the AC category and of the run-length-encoded sequence of zeros
preceding this AC coefficient) of an AC coefficient. In the classic JPEG DNA
experimentation software’s algorithm, the Goldman encoder encodes those
values. To study if and how the proposed variable length coding algorithm
improves the compression, we need to make it so that during the comparison,
the only variant is the choice of the variable length coder for the categories
and run/categories.

To evaluate the performance of the two codecs we computed the compression
ratio, expressed in bits per nucleotide. It corresponds in our case to the
ratio of the size of the input image (in bits) over the total size of the DNA
sequences (i.e., the number of nucleotides) generated to compress the image.
This compression ratio is evaluated for the same quality parameters for the
two different codecs. The PSNR in function of the compression ratio is
presented in the figure 3. The kodak dataset was used to evaluate the
performance of the new image compression algorithm.

Through all the images in the test image dataset, our new variable-length
coder has consistently improved the performance of the compression algorithm
over the one using the Goldman entropy coder. For example with the image
kodim23 from the kodak dataset, the compression rate gain varies between 0.5
and more than 2 bits per nucleotide as shown in the figure 3.

Figure 3: Compression performance comparison for the image kodim23 from the
kodak dataset between the proposed quaternary encoder and the Goldman encoder.
Both encoder have been implemented in the JPEG DNA experimental software, with
$max_{hl}=3$

The relative gain is less important than when using the entropy coder as a
standalone coder because in the case of the image codec, other coders are
used, and those have not been optimized.

##  VI Conclusion

We have developed a new variable-length tree coder adapted to DNA data
storage. This new coder has shown improvements over the Goldman entropy coder
that is usually chosen in variable length DNA coding problems. We have even
shown that this new coder would do better than any ternary entropy coder by
having an expected lenth inferior to the ternary entropy of the source. We
then have integrated this new entropy coder inside an image codec to analyze
its influence on the performance of the codec. We have seen consistent
improvements on the compression rate with the dataset that we used for
testing.

In the future, it would be interesting to study the resilience of this new
encoder to the noise produced by the different processes of the DNA storage
workflow, namely the synthesis, storage and sequencing. Developing new
consensus algorithm, or maybe even error correction codes could also bring a
lot to this new coder. Finally, it could be interesting to mathematically
establish theoretical bounds to the performance of the proposed coder.

## References

  * [1] M. Dimopoulou, E. G. San Antonio, and M. Antonini, “A JPEG-based image coding solution for data storage on DNA,” _EUSIPCO_ , 2021.
  * [2] S. M. H. T. Yazdi, R. Gabrys, and O. Milenkovic, “Portable and error-free DNA-based data storage,” _Scientific Reports_ , vol. 7, no. 1, pp. 1–6, 2017.
  * [3] R. N. Grass, R. Heckel, M. Puddu, D. Paunescu, and W. J. Stark, “Robust chemical preservation of digital information on dna in silica with error-correcting codes,” _Angewandte Chemie International Edition_ , 2015\.
  * [4] Y. Erlich and D. Zielinski, “Capacity-approaching DNA storage,” _bioRxiv_ , p. 074237, 2016.
  * [5] C. Pan, S. M. Hossein Tabatabaei Yazdi, S. Kasra Tabatabaei, A. G. Hernandez, C. Schroeder, and O. Milenkovic, “Image processing in dna,” in _ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2020, pp. 8831–8835.
  * [6] M. Antonini, L. Cruz, E. Da Silva, D. Melpomeni, T. Ebrahimi, S. Foessel, E. Gil, S. Antonio, G. Menegaz, F. Pereira, X. Pic, A. Pinheiro, and M. Raad, “DNA-based Media Storage: State-of-the-Art, Challenges, Use Cases and Requirements,” Oct. 2021, working document of the JPEG DNA AHG. [Online]. Available: https://hal.archives-ouvertes.fr/hal-03447138
  * [7] N. Goldman, P. Bertone, S. Chen, C. Dessimoz, E. M. LeProust, B. Sipos, and E. Birney, “Towards practical, high-capacity, low-maintenance information storage in synthesized DNA,” _Nature_ , vol. 494, no. 7435, p. 77, 2013\.
  * [8] M. Dimopoulou, M. Antonini, P. Barbry, and R. Appuswamy, “A biologically constrained encoding solution for long-term storage of images onto synthetic DNA,” _EUSIPCO_ , 2019.
  * [9] C. E. Shannon, “A mathematical theory of communication,” _Bell System Technical Journal_ , vol. 27, pp. 379–423, 1948.
  * [10] M. Ross, C. Russ, M. Costello, A. Hollinger, N. J. Lennon, R. Hegarty, C. Nusbaum, and D. B. Jaffe, “Characterizing and measuring bias in sequence data,” _Genome Biology_ , vol. 14, no. 5, p. R51.


Certification of MPC-based zonal controller security properties using accuracy-aware machine learning proxies

Pierre HOUDOUIN, Manuel RUIZ, Lucas SALUDJIAN, Patrick PANCIATICI
French Transmission System Operator, RTE, Paris, France

{pierre.houdouin, manuel.ruiz<EMAIL_ADDRESS>
Submitted to the 23nd Power Systems Computation Conference (PSCC 2024).

The fast growth of renewable energies increases the power congestion risk. To address this issue, the French Transmission System Operator (RTE) has developed closed-loop controllers to handle congestion. RTE wishes to estimate the probability that the controllers ensure the equipment's safety to guarantee their proper functioning. The naive approach to estimating this probability relies on simulating many randomly drawn scenarios and then using all the outcomes to build a confidence interval around the probability. Although theory ensures convergence, the computational cost of power system simulations makes such a process intractable.

The present paper aims to propose a faster process using machine-learning-based proxies. The amount of required simulations is significantly reduced thanks to an accuracy-aware proxy built with Multivariate Gaussian Processes. However, using a proxy instead of the simulator adds uncertainty to the outcomes. An adaptation of the Central Limit Theorem is thus proposed to include the uncertainty of the outcomes predicted with the proxy into the confidence interval. As a case study, we designed a simple simulator that was tested on a small network. Results show that the proxy learns to approximate the simulator's answer accurately, allowing a significant time gain for the machine-learning-based process.

Certification of security properties, Congestion management, Multivariate Gaussian processes, NAZA, Proxies

§ INTRODUCTION

Integrating renewable energies on a large scale poses challenges in the operation and management of power systems. It leads to unpredictable and variable flow injections into transmission lines, increasing the risk of power congestion. To address these challenges, the French Transmission System Operator (TSO), RTE, has adopted a decentralized management approach, dividing the entire system into sub-transmission areas (zones). Real-time constraints within each zone are managed through a local closed-loop controller called NAZA [1], [2], [3]. Designed to handle local problems with local actions, it enables the management of battery devices, topological modifications, and curtailment of renewable production [4] inside the zone. Alongside the massive deployment of these controllers across the network, RTE aims to obtain guarantees regarding their proper functioning. Given a set of renewable power production scenarios, RTE seeks to compute the probability of NAZA ensuring the equipment's safety, referred to as $\mathbf{p_{safe}}$ thereafter. Equipment's safety is ensured if congestions are avoided during the scenario simulation. If any line in the zone becomes overloaded during the simulation, it is considered a security threat.

The estimation of $\mathbf{p_{safe}}$ relies on the law of large numbers. The most basic approach is the brute-force process: scenarios are randomly drawn and simulated until $\mathbf{p_{safe}}$ is accurately estimated. Central Limit Theorem (CLT) [5] ensures the convergence of this process. The more iterations are performed, the more accurate the estimation becomes [6]. In our case, only extreme and thus rare-to-observe scenarios will likely pose a security threat. A large number of scenarios must, therefore, be drawn to observe enough threatening situations and obtain a reliable estimation of the probability. Although it theoretically works, the drawback of the brute-force process is that it requires a lot of simulations. As power system simulations are computationally costly, such a process is intractable in this case. However, the outcome of many scenarios can be figured out without a simulation. For instance, a scenario without wind will surely not lead to congestion if the network is correctly exploited. Hence, learning a proxy of the simulator to avoid unnecessary simulations and focus only on the interesting ones can lead to considerable time gain.

The use of proxies to reduce global computational costs meets a growing interest in all domains of power systems. In 2016, Canyasse et al. [7] used supervised learning algorithms to build real-time proxies for solving Alternative Current Optimal Power Flow. Two years later, Duchesne [8] investigated machine-learning proxies to deal with Security Constraint Optimal Power Flow in the context of operation planning. Also, in 2018, Dalal et al. [9], [10] considered the Nearest Neighbors algorithm to predict short-term decision outcomes and applied it to the outage scheduling problem. More recently, in 2022, Chen et al. [11] proposed a deep-learning-based proxy to solve the Security Constraint Economic Dispatch problem and handle real-time applications efficiently.

In this paper, we will also use a proxy to address the computational limitations of the brute-force approach. Our new proxy-based process aims to avoid useless simulations and achieve a faster estimation of $\mathbf{p_{safe}}$. Using batches of simulations, we train a multivariate Gaussian process (MGP) [12] to predict the scenario outcome. Gaussian process (GP) also became very popular in the power systems community [13], [14], [15]. Indeed, they provide exact confidence intervals around the prediction and enable the incorporation of prior knowledge of the system into the proxy [16]. Parameters of the conditional distribution are continuously updated with the new batch of data to enhance its accuracy throughout the process. The confidence interval provided with the GP's prediction enables permanent assessment of the prediction's quality. If the uncertainty in the prediction of the scenario outcome is acceptable, we keep the prediction and avoid a simulation. Otherwise, we perform a simulation. This leads to significantly sped-up iterations when the simulation is not performed. Finally, the CLT is adapted using the Lyapunov version to include the uncertainty of the proxy's predictions in the confidence interval of $\mathbf{p_{safe}}$'s estimation.

The objective of this work is two-fold:

* Train an accurate proxy of a black-box simulator using an MGP

* Show that with an accurate proxy, the machine-learning-based process converges much faster toward $\mathbf{p_{safe}}$ for a given precision, with even a greater gain for demanding precision

The article's structure is as follows: Section II introduces the notations, reminders on the MGP, and problem formulation. Section III presents both the brute-force and the proxy-base processes. Section IV provides the computational results of the case study. Conclusions and future perspectives are discussed in Section V.

§ PRELIMINARIES AND PROBLEM FORMULATION

§.§ Notations

Throughout this paper, upper-case (lower-case) boldface letters will be used for matrices (column vectors), and $(.)^T$ denotes the transposition.
Given a zone of the electric network, we use the following notations :

* $\mathcal{Z}_L=\{L_1,...,L_L\}$ is the set of lines in the considered zone, $L$ is its cardinality

* $\mathcal{Z}_N=\{N_1,...,N_N\}$ is the set of nodes in the considered zone, $N$ is its cardinality

* $P_n$ is the allowed renewable power injection at node $N_n$

* $PA_n$ is the available renewable power injection at node $N_n$

* $P_n^{max}$ is the maximum renewable power that can be produced at node $N_n$, and $\frac{PA_n}{P_n^{max}} \in [0,1]$ is the relative available renewable power injection

* $F_l$ and $\bar{F_l}$ are respectively the power flow and the IST of line $L_l$

§.§ Multivariate gaussian process

Consider an unknown function $h: \mathbb{R}^N \rightarrow \mathbb{R}^L$ that represents the non-deterministic answer $\bs{y} \in \mathbb{R}^L$ of a system to an input $\bs{x} \in \mathbb{R}^N$. We model $\bs{y} = h(\bs{x}) + \mathcal{N}(0, \sigma_0^2 \bs{I})$. Multivariate Gaussian process [17] (MGP) regression aims to learn the underlying dynamic of $h$ by supposing that it is the realization of a multivariate stochastic Gaussian process. The following elements fully characterize the multivariate stochastic process:

* Its mean function : $\mu : \mathbb{R}^N \rightarrow \mathbb{R}^L$

* Its correlation (kernel) function : $k : \mathbb{R}^N \times \mathbb{R}^N \rightarrow \mathbb{R}$

* Its output covariance matrix : $\bs{\Omega}$, a $L \times L$ matrix

Prior knowledge of the system's dynamic can be incorporated into these three elements. For a given $\bs{x} \in \mathbb{R}^N$, $\mu(\bs{x})$ represents the prior expected value of $h(\bs{x})$. Usually, as no specific information about $h$ is known, the choice is $\mu = 0$, which is also our choice here. The kernel function $k(\bs{x_1},\bs{x_2})$ measures to what extent $\bs{y_1}$ and $\bs{y_2}$, the outputs of $\bs{x_1}$ and $\bs{x_2}$, are correlated. Its selection has to be adapted to the problem. For our problem, we choose the stationary squared exponential kernel function :

\begin{align*}
k(\bs{x_1}, \bs{x_2}) = \sigma_f^2 \exp \left( - \frac{|\bs{x_1}-\bs{x_2}|^2}{2l^2} \right) + \sigma_0^2 \delta(\bs{x_1}, \bs{x_2})
\end{align*}

It is a parametric function with parameters $\bs{\theta} = \left[ \sigma_0^2, \sigma_f^2, l \right]$. Parameter $\sigma_0^2$ represents the variance of the system's answer to the same input, $\sigma_f^2$ reflects the signal variance, and $l$ corresponds to the characteristic length-scale of the kernel. The optimal parameters are learned by maximizing the likelihood [12], [18], [19]. The squared exponential kernel results in a smooth prior on $h$ and reflects that similar scenarios are expected to have similar threat potential for the network. Knowing one scenario's outcome only provides local information about neighbors' scenarios' outcomes. The correlation between their output indeed decreases exponentially with the distance. Finally, the process's output covariance matrix reflects the correlations between the output's coordinates. This matrix allows the posterior distribution of the MGP to preserve the covariance structure of the output. The choice of $\bs{\Omega} = I_{L}$ implies that coordinates are decorrelated. The MGP then boils down to training $L$ separate univariate GP to predict each coordinate independently. In this work, $\bs{\Omega}$ is fixed at the beginning of the process. It represents the line's maximum flow correlation and is chosen appropriately regarding the network's structure. An extension of this work to integrate an estimation of $\bs{\Omega}$ across the process's iterations is planned. We now define the following variables :

* $\bs{X_m} = \left[\bs{x_1},...,\bs{x_m}\right]$ the inputs vectors

* $\bs{Y_m} = \left[\bs{y_1}^T,...,\bs{y_m}^T\right]$ the corresponding simulated outputs

* $\bs{M_m} = \left[\mu(\bs{x_1}),...,\mu(\bs{x_m})\right]$ the corresponding mean vectors

* $\bs{\Sigma_m} = \left( k(\bs{x_i}, \bs{x_j}) \right)_{i,j \in [1,m]}$ the covariance matrix between the input vectors

* $\bs{x_{m+1}}$ a new input vector

* $\bs{\Sigma_{[1, m], m+1}} = \left[ k(\bs{x_1}, \bs{x_{m+1}}),...,k(\bs{x_m}, \bs{x_{m+1}}) \right]$ the covariance vector between the new input and the old ones

The joint distribution of $\bs{Y_{m+1}}$ is a matrix-variate gaussian distribution [20] $\mathcal{M} \mathcal{N} \left( \bs{M_{m+1}}, \bs{\Sigma_{m+1}}, \bs{\Omega} \right)$ [17]. Similarly to the univariate GP, we can derive the conditional distribution of $\bs{y_{m+1}} | \bs{y_1},...,\bs{y_m}$ :

\begin{align*}
&\bs{\mu_*} = \mu(\bs{x_{m+1}}) + \bs{\Sigma_{[1, m], m+1}}^T \bs{\Sigma_m}^{-1}(\bs{Y_m}-\bs{M_m}) \\
&\sigma_* = k \left(\bs{x_{m+1}}, \bs{x_{m+1}}\right) - \bs{\Sigma_{[1, m], m+1}}^T \bs{\Sigma_m}^{-1} \bs{\Sigma_{[1, m], m+1}}^T \\
&\bs{y_{m+1}} \in \mathbb{R}^L \sim \mathcal{N} \left( \bs{\mu_*}, \sigma_* \bs{\Omega} \right)
\end{align*}

The obtained multivariate conditional distribution is very similar to the univariate case. The observed points matter for the mean $\bs{\mu_*}$ and the covariance matrix's scale  $\sigma_*$. The covariance matrix's shape is always $\bs{\Omega}$, which ensures that the conditional distribution preserves the correlations between the output coordinates.

§.§ Problem formulation

Given :

* A zone of the power network

* A black-box simulator $\mathcal{S}$ to simulate the zone's behavior in answer to power injections

* A set of renewable power injection scenarios $\mathcal{X}$

, we aim to estimate the probability $\mathbf{p_{safe}}$ to draw a scenario $\bs{x} \in \mathcal{X}$ well-handled by NAZA, i.e., that do not lead to overloaded lines.

Hypothesis on the zone

The zone is composed of $L$ lines, $N$ generator nodes, and a NAZA controller to handle congestion management. We consider the following working hypothesis :

* NAZA controller can curtail each node's renewable energy power injections, thus reducing $P_n \leq PA_n$.

* Each generator produces the maximum available or allowed power

* The loads are constants

* DC approximation: we neglect the reactive power

Renewable energy generation scenario

We model a renewable energy production scenario $\bs{x}$ by the vector of the $N$ available renewable power injections at each network node. To facilitate MGP learning, we draw the relative available renewable power injection. Our set of scenarios $\mathcal{X}$ then boils down to $[0, 1]^N$, which is easier to explore and will be referred to as the information structure. To preserve the spatial correlation of the nodes, the power injection vector $\bs{x}$ is generated using a truncated multivariate Gaussian distribution.

\begin{align*}
\bs{x} = \begin{pmatrix} x_1 \\ . \\ . \\ . \\ x_N \end{pmatrix} \sim \mathcal{N}_{|\mathcal{X}}\left(0, \bs{\Sigma_{nodes}} \right) \longrightarrow \begin{pmatrix} PA_1 = x_1 \times P^{max}_1 \\ . \\ . \\ . \\ PA_N = x_N \times P^{max}_N \end{pmatrix}
\end{align*}

We consider here a very restrictive subspace of low-dimensional scenarios. Our main goal is to show how an MGP can learn the multivariate answer of a black-box function. Future work will deal with more complex scenarios, including a temporal aspect, to simulate the network's behavior on many time steps.

Network simulation

The network's evolution is simulated using a network simulator $\mathcal{S}$. It takes as input a renewable power production scenario $\bs{x}$ and returns the vector of the lines' relative power flow $\bs{y} \in \mathbb{R}^{L}$. For line $l$, the relative power flow is $\bs{y}^{(l)} = \frac{F_l}{\bar{F_l}}$. The network simulator used here is a simple simulator implemented by the authors (more details are provided in the experiments section). Again, an extension of this work, on a real zone of the French transmission network and using a real network simulator [21], is planned.

Well-handled scenarios

A scenario $\bs{x}$ is considered well-handled by NAZA if no congestions on any lines occur. Line $L_l$ is considered congested if its power flow exceeds the IST. Thus, denoting by $f$ the Boolean function that indicates whether a scenario is well-handled or represents a security threat, we have :

\begin{align*}
f\left(\bs{y}\right) = z = \left\{
\begin{array}{ll}
1 & \mbox{if } \forall l \in [1, L], \bs{y}^{(l)} <= 1 \\
0 & \mbox{if } \exists l \in [1, L], \bs{y}^{(l)} > 1
\end{array}
\right.
\end{align*}

[draw] at (0,0) $\begin{pmatrix} x_1 \\ . \\ . \\ . \\ x_N \end{pmatrix}$;
[draw] at (2.7,0) $\begin{pmatrix} P_1^{max} \\ . \\ . \\ . \\ P_N^{max} \end{pmatrix}$;
[draw] at (6,0) $\begin{pmatrix} PA_1 = x_1 \times P^{max}_1 \\ . \\ . \\ . \\ PA_N = x_N \times P^{max}_N \end{pmatrix}$;
[draw,align=center] at (6,-3.2)   Simulator $\mathcal{S}$

Zone $\mathcal{Z}$

Lines features;
[-] (6.8,-3.32)–(7.3,-3.32);
[-] (6.8,-3.32)–(7.1,-3.1);
[-] (7.3,-3.32)–(7.1,-3.1);
(6.8,-3.32) circle[radius=2pt];
(7.3,-3.32) circle[radius=2pt];
(7.1,-3.1)    circle[radius=2pt];
[draw] at (2.6,-3.2) $\begin{pmatrix} y_1 \\ . \\ . \\ . \\ y_L \end{pmatrix}$;
[red, font=] (4.5,-1.6) node Black-box;
[red, thick,dashed] (1.5,-4.7) – (8.1,-4.7) – (8.1,1.5) – (1.5,1.5) – (1.5,-4.7);
[->,ultra thick] (0.65,0)–(1.85,0) node;
[->,ultra thick] (3.55,0)–(4.07,0) node;
[->,ultra thick] (6,-1.18)–(6,-2.13) node;
[->,ultra thick] (4.41,-3.2)–(3.22,-3.2) node;
[draw,align=center] at (0.2,-3.2) $z \in \{0, 1\}$;
[->,ultra thick] (1.98,-3.2)–(1.05,-3.2) node;
(1.7,-3.6) node $f$;

Black-box function taking as input a renewable energy generation scenario $\bs{x}$ and returning $0$ if a congestion occurred, $1$ if security was ensured

The zone, the nodes, lines features, and the simulator constitute a black-box function, as illustrated in Fig.1. The goal of this work is to estimate as precisely and as fast as possible the probability for the black-box function to return $1$: $\mathbf{p_{safe}} = \mathbb{P} \left( f \circ \mathcal{S} (\bs{x}) = 1 \right)$. The following section will detail two processes to estimate $\mathbf{p_{safe}}$.

§ PROCESSES PRESENTATION

§.§ Brute-force process

The brute-force process is the vanilla way to estimate $\mathbf{p_{safe}}$. We randomly sample scenarios in $\mathcal{X}$, compute $f \circ \mathcal{S}$ using the network simulator, and use the CLT to build a confidence interval around $\mathbf{p_{safe}}$. We iterate the sampling until a satisfying confidence interval is reached. Given the $m$ simulation results $z_1,...,z_m$, $\overline{z_m} = \frac{1}{m} \sum_{i=1}^m z_i$ their average, $Q_{\alpha}$ the $1-\alpha$ quantile of the standard normal distribution, the asymptotic confidence interval [6] of level $1-\alpha$ around $\mathbf{p_{safe}}$ is :

\begin{align*}
I_{\alpha} &= [p_{min}, p_{max}]\\
p_{min} &= \frac{\overline{z_m} + \frac{Q_{\frac{\alpha}{2}}^2}{2m} - Q{\frac{\alpha}{2}} \sqrt{ \frac{\overline{z_m}\left(1-\overline{z_m}\right)}{m} + \frac{q_{\frac{\alpha}{2}}^2}{4m^2}}}{1 + \frac{Q_{\frac{\alpha}{2}}^2}{m}} \\
p_{max} &= \frac{\overline{z_m} + \frac{Q_{\frac{\alpha}{2}}^2}{2m} + Q_{\frac{\alpha}{2}} \sqrt{ \frac{\overline{z_m}\left(1-\overline{z_m}\right)}{m} + \frac{q_{\frac{\alpha}{2}}^2}{4m^2}}}{1 + \frac{Q_{\frac{\alpha}{2}}^2}{m}} \\
\end{align*}

The brute-force process is summarized in Algorithm 1.

Brute-force process
Inputs :

* scenario sampler $\mathbb{X}$

* network simulator $\mathcal{S}$

* confidence interval precise enough : $CI_{good}?$ tool

* confidence interval level $1 - \alpha$

Output : $I = [p_{min}, p_{max}]$, $1 - \alpha$ confidence interval level around $\mathbf{p_{safe}}$
[1]
Initialize $m=0$ and $I_m = [0, 1]$
Initialize the bag of outputs $\mathcal{B} = \{\}$
not $CI_{good}?(I_m)$
Draw a scenario $\bs{x_m}$ using $\mathbb{X}$
Increase $m=m+1$
Compute the outcome $z_m = f \circ \mathcal{S}(\bs{x_m})$
Add outcome $z_m$ to the bag $\mathcal{B} = \mathcal{B} \cup \{z_m\}$
Compute $I_m$ of level $1 - \alpha$ using $\mathcal{B}$

§.§ Proxy-based process

The issue with the brute-force process is the very high amount of required simulations to get an accurate estimation of $\mathbf{p_{safe}}$, especially if $\mathbf{p_{safe}}$ is close to 0 or 1. Fig.2 depicts the number of simulations required for the brute-force process to reach the desired relative precision on the estimated probability for different values of $\mathbf{p_{safe}}$. To reduce the number of required simulations, we propose to learn the response $\bs{y} = \mathcal{S}(\bs{x})$ of the simulator to a scenario. We train an MGP-based proxy $\Tilde{\mathcal{S}}$ with the previously simulated scenarios.

precision, Probability_0_5, Probability_0_75, Probability_0_95
35, 59.17, 150.04, 947.69
30, 80.81, 202.22, 1288.96
25, 113.84, 305.92, 1914.76
20, 180.13, 480.24, 3111.04
15, 316.6, 865.35, 5528.19
10, 700.22, 1970.18, 12548.12
05, 2725.44, 7975.07, 50392.28
01, 7102.21, 27910.14, 163123.48

ymode =log,
xlabel=Precision in %,
ylabel=Number of simulations,
xmin=0, xmax=40,
ymin=1, ymax=1000000,
xtick=0, 10, 20, 30, 40,
ytick=0, 10, 100, 1000, 10000, 100000, 1000000,
legend pos=north west,
grid style=dashed,
x dir=reverse

[color=blue, mark=square]
table [x=precision, y=Probability_0_5, col sep=comma] number_of_simulations.csv;

[color=green, mark=square]
table [x=precision, y=Probability_0_75, col sep=comma] number_of_simulations.csv;

[color=red, mark=square]
table [x=precision, y=Probability_0_95, col sep=comma] number_of_simulations.csv;

Number of simulations required for the brute-force process to reach a desired relative precision, depending on $\mathbf{p_{safe}}$ probability, each point is the average of 100 simulations

In this new process, we sample batches of scenarios. For each input scenario $\bs{x_i}$, the MGP-based proxy provides the probability density function (PDF) $g_i$ of the output $\bs{y_i} \in \mathbb{R}^L$. Integrating this PDF allows us to compute the probability $p_i = \int_{[1, \infty]^L} g_i(y) dy$ for scenario $\bs{x_i}$ to lead to congestion. If the confidence in the outcome is high enough, i.e., if $p_i$ is either very small or very high, we keep the prediction and avoid a simulation. Otherwise, a simulation is performed to obtain the outcome. Fig.3 illustrates the decision mechanism to decide when to keep the prediction or simulate the outcome. At the beginning of the process, only a few outcomes are available to train the proxy. Its quality is poor, and the prediction's confidence is often insufficient to avoid the simulation. Most scenarios are thus simulated. However, the more samples are drawn, the better the proxy becomes. At some point, most predictions will be confident enough to avoid the simulation, resulting in a considerable time gain for the probability estimation. The proxy-based process is detailed in Algorithm 2.

[gray,opacity=0.5] (0,0) – (4,0) – (4,1) – (0,1) – (0,0);
[red,opacity=0.5] (0,0) – (0,1) – (-1.9,1) – (-1.9,0) – (0,0);
[green,opacity=0.5] (4,0) – (4,1) – (6,1) – (6,0) – (4,0);
[->,ultra thick] (-1.9,-0.7)–(6,-0.7) node[right]$p_i$;
[-,thick] (-1.9,-0.7)–(-1.9,-0.85) node[below]$0$;
[-,thick] (6,-0.7)–(6,-0.85) node[below]$1$;
[-, green, ultra thick, dashed] (4,-0.7)–(4,1) node[above]$p_{sup}=0.99$;
[-, red, ultra thick, dashed] (0,-0.7)–(0,1) node[above]$p_{inf}=0.01$;
[red, font=] (-1,-0.2) node Prediction is :;
[red, font=] (-1,-0.5) node congestion;
[green, font=] (5,-0.2) node Prediction is : ;
[green, font=] (5,-0.5) node safety ensured;
[gray, font=] (2,-0.2) node Proxy is not trusted : ;
[gray, font=] (2,-0.5) node a simulation is performed;

Proxy use decision mechanism

We compute the confidence interval using only simulated outcomes in the brute-force process. There is no uncertainty on any $z_i$. In the machine learning process, it is different. When the result is predicted, $z_i$ is no longer known with 100% confidence because $0 \leq p_i \leq 1$. Such uncertainty must be considered in the confidence interval around the estimation of $\mathbf{p_{safe}}$. The classic version of the CLT must be adapted to our problem.

Proxy-based process
Inputs :

* scenario sampler $\mathbb{X}$

* network simulator $\mathcal{S}$

* MGP proxy $\Tilde{\mathcal{S}}$

* confidence interval precise enough : $CI_{good}?$ tool

* confidence interval level $1 - \alpha$

Output : $I = [p_{min}, p_{max}]$, $1 - \alpha$ confidence interval level around $\mathbf{p_{safe}}$

Hyperparameters :

* Batch size $m_{batch}$

* MGP parameters

* Proxy usage confidence threshold $p_{inf}$ and $p_{sup}$

[1]
Initialize $m=0$ and $I_m = [0, 1]$
Initialize proxy $\Tilde{\mathcal{S}}$
not $CI_{good}?(I_m)$
Draw $m_{batch}$ scenarios $\bs{x_1},...,\bs{x_{m_{batch}}}$ using $\mathbb{X}$
Increase $m=m+1$
$i \in [1, m_{batch}]$
Compute the probability $p_i = \mathbb{P}(z_i=1) = \Tilde{\mathcal{S}}(\bs{x_i})$
$p_i \leq P_{min}$ or $p_i \geq P_{max}$
Draw $w_i \sim B(p_i)$
Set $q_i = p_i$
Set $q_i = 1 - p_i$
Set $w_i = \mathcal{S}(\bs{x_i})$
Set $q_i = 1$
Add outcome $(w_i, q_i)$ to the bag $\mathcal{B} = \mathcal{B} \cup \{(w_i, q_i)\}$
Update MGP parameters with simulated samples in $\mathcal{B}$
Compute $I_m$ of level $1 - \alpha$ using $\mathcal{B}$

§.§ Central Limit Theorem adaptation

Let $\bs{x_1},...,\bs{x_m} \in \mathcal{X}$ be $m$ scenarios :

* $z_1,...,z_m \in \{0, 1\}$ denote the true answers given by the simulator

* $\forall i \in [1, m], \mathbb{P}(z_i = 1) = \mathbf{p_{safe}}$

* $w_1,...,w_m \in \{0, 1\}$ denote the answers' prediction

* We define $\forall i \in [1, p], q_i = \mathbb{P}(w_i = z_i)$

We suppose that we observe $w_1,...,w_m$ and $q_1,...,q_m$, we look for a confidence interval of level $1 - \alpha$ around $\mathbf{p_{safe}}$.

Proposition 1

We set :

\begin{align*}
\overline{Z_m} &= \frac{\sum_{i=1}^m w_i - (1-q_i)}{\sum_{i=1}^m 2q_i - 1} \\
v_m &= \frac{\left(\sum_{i=1}^m 2q_i-1\right)^2}{\sum_{i=1}^m (2q_i-1)^2} \\
\sigma_m^2 &= \frac{\sum_{i=1}^m q_i(1-q_i)}{\sum_{i=1}^m (2q_i-1)^2} \\
\end{align*}

\begin{align*}
\sqrt{v_m} \frac{\overline{Z_m} - \mathbf{p_{safe}}}{\sqrt{\sigma_m^2 + \mathbf{p_{safe}}(1-\mathbf{p_{safe}})}} \longrightarrow_d \mathcal{N}(0, 1)
\end{align*}

Proof : See Annex A

The MGP does not provide us directly $w_i$, nor $\mathbb{P}(w_i = z_i)$, but only $p_i = \mathbb{P}(z_i=1|MGP \; prediction)$. The previous proposition cannot be directly applied. The issue is overcome by generating $w_i$ using a Bernoulli distribution $B(p_i)$ and choosing $q_i = p_i$ if $w_i = 1$ or $q_i = 1 - p_i$ if $w_i = 0$. Then, we have $\mathbb{P}(w_i = z_i) = q_i$ in all cases. All the information the MGP provides is contained in $q_i$, and the random drawings ensure the correctness of the propositions.

If there is no uncertainty on the observations $w_1,...,w_m$, i.e. $\forall i \in [1, p], q_i=1$, then $w_i \sim B(\mathbf{p_{safe}})$ and we retrieve the CLT for the $z_i$. On the contrary, if predictions provide no information, i.e., $\forall i \in [1, p], q_i=0.5$, then $w_i \sim B(1/2)$ and we retrieve the CLT for balanced Bernoulli variables. Our resulting CLT is a combination of these two CLTs. The balance between the two is determined by the $q_i$. Both terms $\overline{Z_p}$ and $v_p$ can be interpreted as the average $\overline{w}$ and $p$ while $\sigma_p^2$ is an additional variance term. It represents a measure of the quantity of uncertainty contained in the observations. The evolution of this term will be studied in the experiments.

Proposition 2

Denoting by $Q_{\alpha}$ the $1-\alpha$ quantile of the standard normal distribution, the asymptotic confidence interval around $\mathbf{p_{safe}}$ is $[p_{min}, p_{max}]$ with :

\begin{align*}
p_{min} &= \frac{\overline{Z_m} + \frac{Q_{\frac{\alpha}{2}}^2}{2v_m} - Q_{\frac{\alpha}{2}} \sqrt{ \frac{\overline{Z_m}\left(1-\overline{Z_m}\right)}{v_m} +\frac{\sigma_m^2}{v_m} \left[ 1 + \frac{Q_{\frac{\alpha}{2}}^2}{v_m} \right] + \frac{Q_{\frac{\alpha}{2}}^2}{4v_m^2}}}{1 + \frac{Q_{\frac{\alpha}{2}}^2}{v_m}} \\
p_{max} &= \frac{\overline{Z_m} + \frac{Q_{\frac{\alpha}{2}}^2}{2v_m} + Q_{\frac{\alpha}{2}} \sqrt{ \frac{\overline{Z_m}\left(1-\overline{Z_m}\right)}{v_m} + \frac{\sigma_m^2}{v_m} \left[ 1 + \frac{Q_{\frac{\alpha}{2}}^2}{v_m} \right] + \frac{Q_{\frac{\alpha}{2}}^2}{4v_m^2}}}{1 + \frac{Q_{\frac{\alpha}{2}}^2}{v_m}}
\end{align*}

Proof : See Annex B

§ COMPUTATIONAL RESULTS

§.§ Details on the case study

Simulator details

Both processes' performances are tested on a small network of $L=5$ lines and $N=10$ generator nodes. Flows on the lines are computed using the PTDF [25], and $P^{max}$ values, considered as features of the zone. When faraway topological changes occur outside the area, it usually results in modifying the PTDF values. Including topological changes in the scenario $\bs{x}$ would result in an exploding dimension, making the MGP inefficient. Such phenomenon is thus considered a disturbance: at each simulation, we add a small white noise $\epsilon_n$ to the PTDF. Our modeling, however, does not embody significant changes close to the zone, as it would completely redefine the PTDF values. It is one of the limitations of this model.

The actual NAZA controller cannot be integrated into our simulator. Our scenarios have no temporal dimension, and NAZA actions are based on the estimation of the past flow increase. To model fictitious actions, we thus add a $\beta \in [0, 1]$ coefficient that curtails the available renewable power for all nodes.

\begin{align*}
F_l = \beta \sum_{n=1}^N (PTDF_{l, n} + \epsilon_n) x_n
\end{align*}

$\beta$ is selected based on the minimum required curtailment that avoids congestion for scenarios with similar total power injection in the zone. More  specifically, given a set of scenarios $\bs{x_1},...,\bs{x_J}$, we compute the minimum required curtailment that avoids congestion: $\beta_j = \arg\max\limits_{\beta \in [0,1]} \left\{ \beta \sum_{n=1}^N (PTDF_{l, n} + \epsilon_n) x^{j}_n < 1 \right\}$. Then, for another scenario $\bs{x}$, we randomly draw its curtailment value $\beta$ among the set of candidates : $Cand(\bs{x})= \left\{ \beta_j, j \in [1, J], \left| \sum_{n=1}^N x_n - x^{j}_n \right| < \eta \right\}$, with default value to $1$ if the set is empty. Fig.4 illustrates the $\beta$ selection procedure. Again, the goal of this simple fictitious simulator is only to illustrate the proxy-based process. It will be tested with real network simulators, including a real NAZA automaton.

GPA NAZA_curtailment_history
3.45609989580526 0.12911159222701973
1.2728777647215304 0
3.058010578361846 0.02821552155458562
2.522593343372752 0
2.867479625098336 0.012802646169760568
2.716709544585837 0.014927130257798016
2.054045217098228 0
3.7164383189479686 0.09806758590167353
3.5964687754085802 0.06251310352766648
3.361337276654551 0.10186520201607108
0.8517708662169517 0
4.158271538159548 0.1650585598723492
3.6948970854325998 0.014728615795026379
3.80956898558158 0.09030757387798316
4.069920420951497 0.1375656608178203
3.7202628786494967 0.08907873333279348
3.6877921945950956 0.05479642471850443
2.976118104279557 0
3.3711005976288124 0.09641693538291773
3.777024783150401 0.13458634802342284
3.498394234260982 0.08899167213978199
3.3506497022095756 0.0945509386982657
3.3129658981148227 0.010706196057713169
4.010889428614909 0.14357307847321021
3.2937813350043434 0.08570775907986616
3.050792944523695 0.01950331036591102
2.145451557371069 0
3.702673609757794 0.0492920086341081
2.7803154741755693 0.019678499961036322
3.867817747385513 0.08526780083097796
3.413098012931306 0.12461859959804833
3.665883803698083 0.039190601674280456
2.4139769069580304 0
2.0021252430397376 0
3.791288111778803 0.0856911758682366
3.443227566717533 0.11110295955559613
3.4260661244477006 0.04120499440619707
2.942467490614776 0.01421170644649039
3.5780971351466357 0.10605404403379681
3.3509344773259246 0.01970377522752431
2.5368201821542775 0
1.9678631508816538 0
2.2585340314178053 0
4.109866034119287 0.15975013455877185
3.473541770454492 0.06981435731056562
3.357258087301771 0.1386415686419037
4.001093147885151 0.1417499866648151
2.4303273307795616 0
3.7678850975318374 0.05270625651342753
2.2877995161223725 0
3.0472182422081415 0
3.2441313793958337 0.0801431844720234
3.3733071338344525 0.04885389296580622
3.5564561103165384 0.03322770868778879
3.179309896859774 0
3.887809378535794 0.13934630855040098
3.494566942678797 0.012457160961090952
1.825881424246795 0
2.9293400497501914 0.03211637182975868
2.137506457380362 0
3.700412830793491 0.09642797725960228
3.721774134578996 0.07598767901884229
2.924407678576621 0.07235698039079208
3.529036886910988 0.07695289998178678
3.3319623792600126 0.027089943990108667
2.991521046655132 0
1.6198818244682338 0
3.676794893814356 0.04460291891184387
2.209205029823665 0
2.7523694710867908 0
3.395603775159949 0.014466071045003082
3.3219573528968542 0.023264184515382036
3.6531619995424456 0.1271212423663216
3.412901946587933 0.05034466342056887
3.519552633822284 0.032523695907089145
1.7597483246478993 0
2.8148590769499773 0.06276141960778536
3.129597562904957 0.028342908029811653
2.7473753886426473 0
3.085757845566033 0.02215907798113681
3.5462122751538683 0.07312119120921368
2.7086489482817564 0
3.8330916589851993 0.10191930238649649
1.853875225722798 0
3.5373609158716715 0.0675720447152679
2.922649750172052 0
2.392666784025685 0
3.0042460100301946 0.01720901986357748
3.5424221895522523 0.07957225033985282
2.550300805565624 0.023823613624239903
3.172799713076422 0.04301426608470227
3.511275932842013 0.03397008390893552
3.4988216002872874 0.15904457434758945
3.5755891276028002 0.09604791592012396
2.25162568894227 0
3.5256502122047526 0.09348759771967685
2.448776900398565 0
3.3944989196302497 0.09293151398613175
2.8300434068870666 0.04648317700685756
3.427223401199781 0.06959379792157816
2.0762757488543846 0
3.181471201270881 0.01000851723908025
3.9729183060188453 0.16889228096500952
3.536610028377611 0.1319219789344408
1.7442002057167996 0
2.8750329152696477 0.015084662171860441
3.2485066057026133 0.06203032992979629
3.6731947631007404 0.09012188840672075
2.137810550546782 0
3.599380032912273 0.0862289415135451
3.1129136459070086 0.022975982946049633
2.629393279068993 0.0340101956584179
3.43641304292506 0.02344634493766351
2.4864793856175864 0
3.47168379113127 0.020483231613795838
2.1794474238636035 0
3.986078666906882 0.07063807802768118
1.7799138057540649 0
2.4189725534300637 0
3.9720201529224557 0.1349978862619644
3.7976254241702274 0.05383286986428093
3.0891359462621617 0.07881215467400027
3.320501253377971 0.013968258887490204
3.5621628598556025 0.07471891460585267
2.7679428631484426 0
2.4777984523517445 0
3.292815080064316 0.024248597967649443
3.602126405192627 0.14369891371730203
3.234855649407157 0.06663513861923198
3.638119908646951 0.02252109131758384
3.78942906031918 0.02256679490452218
2.554511786742309 0
3.919000299247643 0.13384621382872552
3.5937306071303134 0.10843177543148574
4.087728296279282 0.11749708246112534
4.318382147892403 0.16678422097587553
3.6449856014477904 0.10167987219998609
4.498842229828157 0.21415870697542017
3.5695718525329756 0.04814042632965421
3.3085188736129556 0.09712121078569402
3.093205981682012 0.025255190716305997
1.7559768383725343 0
3.5438481210298765 0.02882867261400024
3.141886065180943 0.05539676515839176
3.66237576612406 0.024965171026552047
4.297142744471975 0.14060370318425708
3.170897664626502 0
2.604344303358398 0
2.907894309936549 0
2.6621645003015426 0
3.1236020164225113 0.03743750611509822
3.7699275835547605 0.055721103510916636
2.4394200224970697 0
2.656693596889232 0.031575436635923104
3.2239821777185833 0.01009005415157771
3.6266146632512934 0.08884533585519383
3.510014923849365 0.0589529479212636
3.5949494853483115 0.06549339115629871
3.0660615149110777 0.03230778155033842
4.112346783075827 0.1288001489440772
3.9635134440911135 0.11217388842322285
3.2240517877459585 0.03245044616361226
2.8074206329123497 0
0.9772047176330906 0
2.312987845129301 0
3.147026522792637 0.07436933233798615
2.65471779218983 0
2.7164961488520643 0
3.695980255183787 0.12619499594980288
3.1752892822379613 0.03249907003403152
3.060559153153824 0.023464380990127043
3.422390100635617 0.06329072495085697
2.9533257211203328 0.062300473839038695
3.7997914538563435 0.06891039214846373
2.2674077301729105 0
3.7371048510814386 0.0749452843024846
3.2460199742367597 0
3.1666289249981787 0.10505820912026453
3.5445452442398167 0.07704892819195275
3.1629974708035844 0.06210181190491959
2.508299273660618 0
1.7270684923607766 0
3.2461962847379118 0.013438650725677515
3.367986812716227 0.025515712649700206
1.8887594858753192 0
3.1286586690902825 0.02961222211218628
3.1794151853890624 0.02637171796609261
1.9588546031823646 0
3.5655158840142804 0.03330724292420706
3.11105506525469 0
3.16107066521483 0.05684463327564959
3.63134138006423 0.05936138413925185
2.037810014678732 0
3.639736440932845 0.01768038500915936
4.192427160098913 0.1924387918066628
3.2192009827974295 0.0656702845259729
3.4002043834337 0.0981105047565427
3.7727487813387683 0.14701744666659444
3.398049349862884 0.017774510742452954
3.069582869207732 0.01649971014870899
1.6538330886030552 0
3.758880151058979 0.02707007215284918
3.0039218379257315 0.03415459286400768
2.160646721180698 0
3.099414357638617 0.05663379481000608
3.288766814954862 0.017424133989517254
1.9842581001254607 0
2.796549913099456 0.050404312553405384
3.022776197521014 0.028729151569844236
3.290288851014314 0.019747748117464492
3.495699480269199 0.06789401588868782
3.2043709183914233 0.04588981944297088
2.027848054318182 0
3.001063574340569 0.024083555580367455
3.2976665709373236 0.04913673413546138
3.4656945665351566 0.03439458952507457
4.129467963967726 0.1328848697479026
3.9565450640631337 0.08990993715193496
3.637490890758924 0.07984425128366224
3.355401363255972 0.017267171321340258
3.0036270789165815 0.05601032777156789
2.8833172731173375 0
1.2818384654751065 0
3.233112552415446 0.05707648450805858
3.3726682125350114 0.019324310874135442
3.482781817909027 0
3.299175437093732 0.13203143942790774
2.070356731343196 0
2.9865326528419236 0.015849859114960507
3.5368271864491385 0.05278127721807302
0.9215811770453575 0
3.2863726218728306 0.05994528514393238
3.4857803853668576 0.07659078234461114
2.970997344054952 0.012300237343501497
2.531240376722203 0
3.7077118761611 0.14247872014969654
3.049499473172542 0.012045864855243282
3.8621873644459477 0.12918092966989547
3.0990659952560415 0.051088619867434626
3.7859312198476482 0.15314209062354545
3.3149721044169898 0.020519812327698728
1.6797423857944556 0
3.002678578893356 0
3.7554110568910195 0.11592883209400695
2.562525471497732 0
3.2700497589532667 0.01640112624711343
2.985425546098975 0
3.066489296978508 0.026625696180899605
3.6407701617554524 0.11616375899089815
1.8713274137948959 0
3.446246300807326 0.061364950672579965
3.9714186938241935 0.12939966134538072
3.7573265296623046 0.09551788560369923
2.869935654908709 0.02677728515391742
1.6109153704150947 0
1.8154493770598852 0
3.073426372111371 0.02681302690163956
2.2488556048396195 0
1.7270948187746844 0
4.06240038774378 0.20449740695053797
3.157239845963983 0.0234272220418269
3.8337660918823837 0.14432065562195184
2.966033653615959 0.026694005470842663
2.9880688855036426 0.012037072693064421
3.6803059446793487 0.09554199988584744
2.861868037429442 0.0645439789888177
3.96895279494376 0.12127211528230841
2.4997981126510713 0
3.410901253379081 0
3.421702905194302 0.009400908272958364
4.314317373548554 0.19023416852089725
1.9909162454697147 0
2.709987178457138 0
2.1957200847957075 0
3.1787518862306507 0.02450251100761025
1.0636879043516612 0
3.337368975343292 0.05176232820507698
3.7064784337720282 0.08694139623163277
3.5289632826923216 0.06615644082490502
3.823750444286141 0.10250542384009309
3.7474580941579054 0.11473242562446928
3.5327833533638957 0
3.6160101350572598 0.08429704456224862
2.395594344731742 0
3.9352004374398146 0.13940343514330067
3.953649829763112 0.083502277354494
3.194135120761309 0.07703560937071974
2.4431511564858717 0
3.115592332419094 0
3.550861743021364 0.02657234787368512
3.4328127594224545 0.06089648800878733
3.5832928168799785 0.09117587010090955
3.3034140670553027 0.028535703418244962
2.5005450905756987 0
3.472408699347136 0.06056477997495724
2.9509385432646558 0
2.4259006376412406 0
2.1198073443428638 0
3.551313347418404 0.13045491537059095
2.222098154377162 0
2.6910205989716127 0
2.2143607442083812 0
2.333529365714691 0
3.150794016345674 0.06359817823427827
1.7355838601716904 0
3.523926167038258 0
2.175670686149708 0
3.097128478839693 0.08548363220367561
3.13817748603517 0.07545775972035285
3.3771733828760944 0.028368725644649984
1.9313415635761417 0
2.4252957517459377 0
3.4921484965207723 0.03815629449770108
2.1207162357661775 0
2.0552202672255815 0
3.943090453859966 0.1451150305548281
3.607778206845919 0.12455156413662993
3.556191941277058 0.04271344217808504
2.4889362349588025 0
4.1720031016190555 0.1348338099909978
3.4740886400086817 0.09344552123686323
3.727436597675051 0.07719477264523666
2.2099572712460325 0
3.460348429324778 0.01462759269827707
3.5847295163579225 0.06031881845965459
3.4558680592279964 0.09667483399377297
3.492442634027688 0.07625648365487647
3.473972631091954 0.08280944126882639
2.812456709201745 0
4.351357712444712 0.17955772872630915
3.917710657172815 0.14601438459042215
3.5277726744739235 0.026435187724507703
3.526759908300651 0.06542002941093983
3.205200235478343 0.053326904029556976
2.888624632326882 0
2.614466700119754 0
3.6731003276341907 0.016229152163297678
3.812229867597117 0.04695069244308514
3.9318112895378623 0.06601019740338676
1.9324065281677205 0
2.424469768623418 0
2.6659644114848846 0.029077568269759713
3.2950196572751285 0.09994324874121996
3.4121959338044214 0.07249580223955887
3.0606029738413616 0.007005573389357811
4.394904441888303 0.169845184750697
2.176705074504522 0
3.7351570859192478 0.0161385350998845
3.9869277331898045 0.1249051418217435
2.0288739749762112 0
2.611423759381487 0
3.5149413925112487 0.04079915359036246
3.1790083179744757 0
3.5205069651102217 0.018339594243795054
2.4564089045765076 0
3.310122159318285 0.05235754413925042
3.3057722884909944 0.010450208686906493
3.722049157623925 0.08528138339708369
3.351599209353997 0
3.6367857647405466 0.11977168717592557
3.103366118761225 0.025816311026043537
3.6978925524771933 0.07708608591078714
2.013216265248665 0
3.7133454609244696 0.07543185357020576
2.802266708077526 0.022876807954357425
3.730218763555091 0.11326625531596746
3.220284920795092 0
3.5302117647964115 0.05913655727002243
3.7536350284680275 0.11438752750409187
4.154216657297665 0.12300861850765621
2.459385318495327 0
3.853461224908841 0.13166126322722094
3.286096976899041 0.1419781686987207
4.06398019015815 0.14649966761137667
3.125696671396525 0.024437577449039183
3.7348184970759624 0.12008932912072967
2.8981983393085367 0.05773496289398461
2.900301063742096 0
2.2010002480209936 0
3.3793372568930686 0.07840576092244611
3.2767688538163022 0.022888246348801156
3.191254570330856 0
1.6380445637102512 0
2.4913858008370964 0
3.768543739333323 0.10110235085890784
3.3370578346076343 0
3.191917081796605 0.023170239143428417
3.9189279780236577 0.11513240230281019
3.6419564029477716 0.08761060550587385
3.7624014409633038 0.07425373485716069
3.2669376408636595 0.04950536641656922
2.0492980836595263 0
2.8647154333641973 0.03951308227460326
3.3395657429139227 0.02758066215198618
2.606733380513811 0
3.5608991990878964 0.0710170769351619
4.219549927927105 0.15372679210108664
1.7509661117655728 0
2.440056926515741 0
3.3923730464211848 0.025767958899749366
3.2766870980551523 0.012847031636705974
3.593674487037431 0.019731073662352938
2.779131506820737 0
3.8710027453425466 0.10563483449027966
4.26430804909753 0.14473916322458724
3.1033546174286464 0.05265091979557879
3.196113510077069 0.050522218071253566
2.8160139174474237 0.029367571850104635
4.176970905774806 0.14726589726889716
2.14987597470375 0
2.852152291797342 0
4.107661508880506 0.19372694279960523
3.357233240991246 0.04846488852034503
1.760229103654868 0
3.2782555819237187 0.055530720880170525
0.7876467012146916 0
3.820461325229816 0.0475868881544903
3.6413493503655596 0.1551154486177545
3.648173969185608 0.0788820165051719
3.8038123205858496 0.07685903598087684
2.6094620443675725 0
2.9835910543186484 0
2.6170858851403778 0
3.34966849265423 0.019653644016832228
3.578147661414905 0.12368421000726026
2.8086843180161774 0.04507362863640296
3.1295310223765296 0
3.9837895555140888 0.127410437360586
1.937434353441033 0
3.228944667759849 0.05047709711846726
3.184573770986045 0.05459792104124106
2.7819022211335054 0
3.5663094222582297 0.024521933449116795
3.2580233385379436 0.029202156133897095
1.9680171494132483 0
3.6796018670819484 0.006226672252120147
2.046685208642858 0
3.537399277759579 0.06586012642109873
3.064206702206582 0.02085256320046258
1.8028357450235593 0
3.706793264368974 0.08950261927131431
3.3518035737283722 0.10546376846224159
3.0427496627219615 0.022889843117741226
2.766389435574233 0.03150925228935797
2.9231759420017096 0.013165229352927965
3.5402423689203593 0.07152836915819771
2.9951477642229993 0.024793505544591588
2.183832536224003 0
2.8337403903006066 0
3.5259605999372834 0.07298010465635019
1.7124771258149454 0
2.6964342162066326 0.02311730310379323
3.163245983456502 0
3.631373295973208 0.04544324683302248
2.740548541330104 0.013103795374666438
3.413178310925967 0.052170490489995336
3.82932271567414 0.10832391844614986
3.1185089267949526 0
4.18831067035903 0.14968318567065542
2.457517521959254 0
1.927659702057512 0
2.3053502774249397 0
3.446801856223955 0.05150286296678297
2.271153172687246 0
3.3520959030738484 0.025409955940751017
2.6450856186319887 0
4.3644004844423 0.1580117767321687
4.176288038758477 0.13790176369292928
2.6762128802946807 0
2.0670010755466546 0
3.7285284257912408 0.12645289611689492
3.7391544409312494 0.06617704018535042
3.4162480408788967 0.0550842749820285
2.9905871032899243 0.021258675949614392
3.073345759069258 0.08092914547064001
2.4017100986231883 0
2.1353758005990953 0
2.558118831296923 0
1.5849565516035022 0
3.175469258623501 0.042107649566042515
2.17283213318528 0
3.12896419441692 0
2.2219849271216123 0
3.649469656676277 0.11042128209707358
2.349377276128018 0
3.7751128102094578 0.04469706590554423
3.301370075105016 0.10904727006647019
3.7129990221325064 0.0807870897578381
3.417300962601692 0.09619474101259962
3.6541163012162383 0.08483223041892418
3.1792296587811744 0.06973356232889705
2.3776447701352645 0
2.7706724947558845 0.04998929340476144
3.972247358416882 0.14224319078317182
3.4993078022169763 0.04026307975262389
2.674848099325873 0
1.7065426957621486 0
3.5341370037720043 0.06969073817901406
3.99263517788724 0.17329147533325473
3.32580748995994 0.0318228354212734
1.9842797039771591 0
3.316812281431669 0.045019266913820855
3.6651103147815123 0.10721938873237491
3.544555353780063 0.06723392601040212
3.603284004483425 0.05129506321604749
3.5365062012727773 0.16887070452447672
2.440836334487818 0
3.8967012581779397 0.18044716348267753
2.738778659150327 0
2.444546208280675 0
3.2331180278128366 0.039054721858027164
1.938340987054424 0
2.6720892998438024 0
3.818927809140456 0.11796816543454969
3.47366467888648 0.022322426002375217
3.0970050462274132 0.015477317321648032
2.1430988503462522 0
2.1862836349882953 0
3.6966618485195433 0.10022553023757241
3.6576846627174 0.0764606285289623
3.2116341482112234 0.03713652217025609
3.060485414624125 0.056309131552384484
3.844918790069135 0.15645673178112873
3.4310690247383957 0.028022022480854725
3.755181641536984 0.12046677257922833
2.126690322171344 0
2.92045252047031 0.019780669546237806
3.1344107241723336 0.013239766884855841
3.176964731124005 0.015811952684519848

xlabel=$\sum_{n=1}^N \bs{x^{(n)}} = $ Total relative injection,
ylabel=curtailment value $\beta$,
xmin=0, xmax=5,
ymin=0, ymax=0.2,
xtick=0, 1, 2, 3, 4, 5,
ytick=0, 0.04, 0.08, 0.12, 0.16, 0.20,
legend pos=north west,
grid style=dashed,
+[only marks, scatter, mark=o, mark size=1.2pt]

[thick, red] (3.8,0) – (3.75,5.7);
[draw,align=center, red] at (2.5, 4.5) $\sum_{n=1}^N \bs{x^{(n)}}$;
[thick, dashed, green] (3.5,0) – (3.5,5.7);
[thick, dashed, green] (4,0) – (4,5.7);
[<->][thick] (3.5,4.5)–(4,4.5) node;
[->][thick] (2,3)–(3.25,2.5) node;
[draw,align=center] at (1.27,3) Random
[black, font=] (3.75,4.75) node $2 \eta$;
[green] (2.75,3.75) node $Cand(\bs{x})$;

Selection of the curtailment value $\beta$ for a given scenario based on historical curtailment decisions

Convergence criterion

We sample random scenarios until a certain number of simulations have been reached. For the brute-force process, the number of simulations matches the number of sampled scenarios. More iterations can be performed for the proxy-based process, depending on the proxy's accuracy.

Computational details
Over the iterations, we might have tens of thousands of simulated scenarios to consider when training the MGP. It, however, implies inverting a covariance matrix with a dimension equal to the number of simulated scenarios, thus over tens of thousands. It is too time- and space-consuming to be efficiently computed. Nevertheless, when computing the likelihood of a new scenario, many previously simulated scenarios will be very different. Their correlation with the new input will be negligible. All these scenarios could thus be neglected when estimating the parameters of the conditional distribution. Inspired by [22], we consider only a small proportion of the simulated scenarios when computing the likelihood of a new input $\bs{x_0}$.

* At most $N^*$ scenarios

* Only neighbors scenarios that are at a distance less than $d^*$: $d(\bs{x}, \bs{x_0}) < d^*$

The last point that deserves to be highlighted is the computation of the integral for $p_i$. Both processes are implemented in Python; we use the Scipy library [23] to approximate the integral numerically. In small dimensions, the evaluation of the integral is fast. In higher dimensions, the time quickly increases, cf Fig.5. It is another limitation to overcome to generalize this approach to bigger networks.

dimension, CDF_computation_time
5, 0.00523949294090271
10, 0.022566256999969482
15, 0.042873699426651
20, 0.05942802143096924
25, 0.07479980988502502
30, 0.09069558310508727
35, 0.09590010261535645
40, 0.11744982137680053
45, 0.14116370482444765
50, 0.15945626530647278
55, 0.17481822619438173
60, 0.1953198261260986
65, 0.21154743556976321
70, 0.23230545954704285
75, 0.25365730319023133
80, 0.2750231835842133
85, 0.29144797101020814
90, 0.29297260951995846
95, 0.29915142908096313

ylabel=Time in $s$,
xmin=0, xmax=100,
ymin=0, ymax=0.3,
xtick=0, 20, 40, 60, 80, 100,
ytick=0, 0.05, 0,1, 0.15, 0.2, 0.25, 0.3,
legend pos=north west,
grid style=dashed,

[color=black, mark=square]
table [x=dimension, y=CDF_computation_time, col sep=comma] CDF_computation_time.csv;
Evaluation time

Integral evaluation time depending on the dimension, 10000 averaged simulations per point

§.§ Proxy performances

In this subsection, we illustrate how the MGP manages to learn the behavior of the black-box function to act as a proxy of the simulator. We plot three different indicators.

First, we display in Fig.6 the evolution of $\sigma_*$ across the process's iterations. $\sigma_*$ is the scale of the covariance matrix of the posterior distribution. A small  $\sigma_*$ reflects a low variance and thus a high confidence in the predicted outcome. $\sigma_*$ entirely depends on the point we wish to predict, rendering the plot irregular. To catch the global trend, we display a moving average.

iteration, sigma_star
1, 1.0
5000, 0.6139041438941205
10000, 0.42650354170738103
15000, 0.33864372737987425
20000, 0.2850107469849657
25000, 0.21679900907328398
30000, 0.18529028647938317
35000, 0.18963324207263593
40000, 0.1410464811332522
45000, 0.14470326208622253
50000, 0.12048875737958457
55000, 0.10696850246791678
60000, 0.12021302153554153
65000, 0.1062675261332601
70000, 0.10030304189009089
75000, 0.09260946312694854
80000, 0.08481394281717529
85000, 0.08192535387564574
90000, 0.06950026190538866
95000, 0.07374680560743685
100000, 0.06699225953998941
105000, 0.06256583362212935
110000, 0.057613079605236275
115000, 0.06223027675917532
120000, 0.062143573072945255

xmin=0, xmax=125000,
ymin=0, ymax=1,
xtick=0, 20000, 40000, 60000, 80000, 100000, 120000,
ytick=0.2, 0.4, 0.6, 0.8, 1,
legend pos=north west,
grid style=dashed,

[color=blue, mark=square]
table [x=iteration, y=sigma_star, col sep=comma] sigma_star.csv;

Evolution of the 200 mobile average of $\sigma_*$ during the proxy-based process

In the beginning, the MGP has no observations from which to learn. The variance of the posterior distribution is $\sigma_*=\sigma_f + \sigma_0=1$, an arbitrary prior choice. Throughout the process, we simulate more scenarios, increasing the chances of having many correlated points close to the new sampled observation. Mathematically, this boils down to a higher value for the Mahalanobis term $\bs{\Sigma_{[1, m], m+1}}^T \bs{\Sigma_m}^{-1} \bs{\Sigma_{[1,m], m+1}}^T$, and thus a lower variance. The decrease of $\sigma_*$ illustrates the increasing density of simulated scenarios in the information structure, directly linked to a downward variance posterior distribution and a high confidence in the prediction. $\sigma_*$ could amount to a loss function, whose decrease attests that our MGP indeed learns the simulator's answer.

The second indicator is related to the evolution of confidence in the outcome's prediction across the process's iterations. We plot in Fig.7 the evolution of the prediction's entropy: $H(z_i) = \frac{-p_i log(p_i) - (1-p_i) log(1-p_i)}{log(2)}$. A small entropy characterizes a probability close to $0$ or $1$, thus a confident prediction, while a higher entropy stands for a more uncertain prediction. Again, we display a moving average to avoid a very irregular plot.

iteration, entropy
1, 1
5000, 0.886580602122271
10000, 0.7818007323666759
15000, 0.7010086986025696
20000, 0.653508115513733
25000, 0.6328702407023237
30000, 0.598842455422774
35000, 0.5463405047599799
40000, 0.537857377473457
45000, 0.49327310756599974
50000, 0.4610095659754656
55000, 0.47238428063042884
60000, 0.48449263336566184
65000, 0.4276415514442126
70000, 0.45579043526739776
75000, 0.3916609189383965
80000, 0.38413670246557174
85000, 0.36293072930423234
90000, 0.34016473840741296
95000, 0.36598621499190037
100000, 0.3669173320927011
105000, 0.33419464048123104
110000, 0.3344703594692345
115000, 0.30412496846895826
120000, 0.3127838728025358

xmin=0, xmax=125000,
ymin=0, ymax=1,
xtick=0, 20000, 40000, 60000, 80000, 100000, 120000,
ytick=0.2, 0.4, 0.6, 0.8, 1,
legend pos=north west,
grid style=dashed,

[color=blue, mark=square]
table [x=iteration, y=entropy, col sep=comma] entropy.csv;

Evolution of the 200 mobile average of the proxy's prediction entropy during the proxy-based process

In the beginning, the posterior distribution is $\mathcal{N}(0, \bs{\Omega})$. $\bs{\Omega}$ scale is adapted so that the prior probability is unbiased: $p_1=q_1=\frac{1}{2}$ and $H(z_1)=1$. The congestion probability is mainly influenced by whether or not $\mu_*$ is close, relatively to $\sigma_*$, to the congestion limit, i.e., the boundaries of the $[0,1]^N$ volume. As $\sigma_*$ decreases, it becomes less and less likely to be close to the border. Thus, the computed probabilities are increasingly extreme, close to $0$ (confidence in an overload) or very close to $1$ (confidence in a safe line).

Finally, the last interesting factor to display is the % of simulated scenarios since the beginning of the certification process. The lower the % is, the more simulations were avoided, and thus the better the proxy is.

iteration, perc_of_simulations
1, 1.0
5000, 0.9996000799840032
10000, 0.9971002899710029
15000, 0.9938670755282981
20000, 0.9887505624718764
25000, 0.9826806927722891
30000, 0.9744675177494083
35000, 0.9657438358904031
40000, 0.9558761030974225
45000, 0.944779004911002
50000, 0.9353812923741525
55000, 0.9256013527026782
60000, 0.9158180696988384
65000, 0.9047553114567468
70000, 0.8944729361009132
75000, 0.8846682044239415
80000, 0.875114061074237
85000, 0.8653780543758313
90000, 0.8563460405995492
95000, 0.8473279228639702
100000, 0.8379216207837924
105000, 0.829411148465253
110000, 0.820628903373606
115000, 0.812192937452718
120000, 0.8042099649169593
125000, 0.7965936272509823

ylabel=% of simulations,
xmin=0, xmax=125000,
ymin=0.75, ymax=1,
xtick=0, 20000, 40000, 60000, 80000, 100000, 120000,
ytick=0.75, 0.8, 0.85, 0.9, 0.95, 1,
legend pos=north west,
grid style=dashed,

[color=blue, mark=square]
table [x=iteration, y=perc_of_simulations, col sep=comma] perc_of_simulations.csv;

Evolution of the total % of simulated scenarios since the beginning during the proxy-based process

Initially, the proxy's confidence is insufficient to replace an accurate simulation, and most scenarios are simulated. Then, as the proxy quality improves and the prediction's entropy decreases, the proxy is used more and more often. Over the last $1000$ of sampled scenarios, less than $50\%$ are simulated, leading to $20\%$ of avoided simulations.

§.§ Processes performances

To conclude the experiments section, we present the performances of each process for probability estimation. Given a maximum amount of simulations, we compare in Table 1 the estimated $\mathbf{p_{safe}}$ probability and the confidence interval obtained by both processes.

The brute-force process performs 100,000 iterations and reaches a relative error of 0.05%. Comparatively, the proxy-based process gets over 125,000 iterations and avoids more than 25,000 simulations. The proxy-based confidence interval is computed with more iterations, leading to a smaller relative error and a smaller length. The performance gap is almost $20\%$ above the brute-force process. Since the proxy is continuously enhanced across the simulations, the gap between the two processes would be even more significant for a higher maximum amount of simulations.

3|c|blue!20 $\mathbf{p_{safe}} = 0.9044$

black     red!25 Brute-force     green!25 Proxy-based

Number of iterations     red!15 125575     green!15 100000

$p_{min}$     red!15 0.9018     green!15 0.9023

$p_{max}$     red!15 0.9067     green!15 0.9064

Relative error     red!15 0.05%     green!15 0.041%

Confidence interval length     red!15 0.0049     green!15 0.0041

Processes performances comparison

[->, thick] (0,0)–(8,0);
[-, green, ultra thick] (1.8,-0.25)–(1.8,0.25);
[-, green, ultra thick] (6,-0.25)–(6,0.25);
[-, green, ultra thick] (1.8,0)–(6,0);
[black] (4,0.5) node $\mathbf{p_{safe}}=0.9044$;
[green] (1.2,-0.6) node $p_{min}=0.9023$;
[green] (6.6,-0.6) node $p_{max}=0.9064$;
[-, red, ultra thick] (1.3,-0.25)–(1.3,0.25);
[-, red, ultra thick] (6.4,-0.25)–(6.4,0.25);
[-, red, ultra thick] (1.3,0)–(1.8,0);
[-, red, ultra thick] (6,0)–(6.4,0);
[-, ultra thick] (4,-0.25)–(4,0.25);
[red] (1.1,0.5) node $p_{min}=0.9018$;
[red] (6.8,0.5) node $p_{max}=0.9067$;
[-, red, ultra thick] (1,1.5)–(1.6,1.5);
[-, green, ultra thick] (1,2)–(1.6,2);
[black] (3.3,1.5) node : Brute-force process;
[black] (3.375,2) node : Proxy-based process;

Confidence intervals obtained by both processes

§ CONCLUSION

In this work, we proposed two processes to estimate the probability of observing congestion on a network. The network and its simulator are a black-box function that outputs the maximum flow observed on each line for a renewable energy production scenario. The brute-force process draws and simulates scenarios until it reaches a satisfying probability estimation. As power system simulations are often very long, they are time-consuming. The proxy-based process trains a Multivariate Gaussian Process to predict the black box's answer, with an exact confidence interval on the prediction. The simulation is avoided when the proxy is confident enough to make the prediction. For a given budget of simulations, it leads to a much more accurate probability estimation, with an error reduced by $20\%$. Many simulations are not performed, allowing the process to run more iterations. The MGP manages here to efficiently learn the multivariate answer of a black-box function with correlated components while providing explicability and theoretical guarantees about the prediction's error. Multivariate Gaussian processes are of great interest for approximating power systems behaviors and can be applied to many problems.

We supposed in this work that the output's correlation matrix $\Omega$ was known, while in practice, it has to be learned. A first extension of this work will be to propose an estimation method of this matrix during the process. Also, the kernel choice can be discussed. This paper aimed to illustrate the proxy-based certification process on a simple black box involving low-dimensional inputs. We aim to extend it to more realistic simulations with a real network and complex scenarios.

[1]
Hoang, Duc-Trung and Olaru, Sorin and Iovine, Alessio and Maeght, Jean and Panciatici, Patrick and Ruiz, Manuel, Power Congestion Management of a sub-Transmission Area Power Network using Partial Renewable Power Curtailment via MPC, 2021 60th IEEE Conference on Decision and Control (CDC)

[2]
Hoang, Duc-Trung and Olaru, Sorin and Iovine, Alessio and Maeght, Jean and Panciatici, Patrick and Ruiz, Manuel, Predictive Control for Zonal Congestion Management of a Transmission Network, 2021 29th Mediterranean Conference on Control and Automation (MED)

[3]
Straub, Clementine and Olaru, Sorin and Maeght, Jean and Panciatici, Patrick, Zonal Congestion Management Mixing Large Battery Storage Systems and Generation Curtailment, 2018 IEEE Conference on Control Technology and Applications (CCTA)

[4]
Meyer, Bruno and Astic, Jean-Yves and Meyer, Pierre and Sardou, François-Xavier and Poumarede, Christian and Couturier, Nicolas and Fontaine, Mathieu and Lemaitre, Christian and Maeght, Jean and Straub, Clémentine, Power Transmission Technologies and Solutions: The Latest Advances at RTE, the French Transmission System Operator, IEEE Power and Energy Magazine

[5]
Fischer, Hans, A history of the central limit theorem: from classical to modern probability theory, 2011 Springer

[6]
Jeffrey T. Morisette and Siamak Khorram, Exact Binomial Confidence Interval for Proportions, 1998 PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING

[7]
Canyasse, Raphaël and Dalal, Gal and Mannor, Shie, Supervised learning for optimal power flow as a real-time proxy 2017 IEEE Power & Energy Society Innovative Smart Grid Technologies Conference (ISGT)

[8]
Duchesne, Laurine and Karangelos, Efthymios and Wehenkel, Louis, Using Machine Learning to Enable Probabilistic Reliability Assessment in Operation Planning, 2018 Power Systems Computation Conference (PSCC)

[9]
Dalal, Gal and Gilboa, Elad and Mannor, Shie and Wehenkel, Louis, Unit Commitment Using Nearest Neighbor as a Short-Term Proxy 2018 Power Systems Computation Conference (PSCC)

[10]
Dalal, Gal and Gilboa, Elad and Mannor, Shie and Wehenkel, Louis, Chance-Constrained Outage Scheduling Using a Machine Learning Proxy 2019 IEEE Transactions on Power Systems

[11]
Wenbo Chen and Seonho Park and Mathieu Tanneau and Pascal Van Hentenryck, Learning optimization proxies for large-scale Security-Constrained Economic Dispatch, 2022 Electric Power Systems Research

[12]
Rasmussen, Bousquet, O., von Luxburg, U., Rätsch, G,  Gaussian Processes in Machine Learning, Advanced Lectures on Machine Learning. ML 2003. Lecture Notes in Computer Science, vol 3176. Springer, Berlin, Heidelberg

[13]
Kocijan, Jus, Modelling and Control of Dynamic Systems Using Gaussian Process Models, 2016 Springer

[14]
Jain, Achin and Nghiem, Truong and Morari, Manfred and Mangharam, Rahul, Learning and Control Using Gaussian Processes 2018 ACM/IEEE 9th International Conference on Cyber-Physical Systems (ICCPS)

[15]
Yandong Yang and Shufang Li and Wenqi Li and Meijun Qu, Power load probability density forecasting using Gaussian process quantile regression 2013 Applied Energy

[16]
Ana M. Ospina and Yue Chen and Andrey Bernstein and Emiliano Dall’Anese,  Learning-based demand response in grid-interactive buildings via Gaussian Processes, 2011 Electric Power Systems Research

[17]
Zexun Chen and Bo Wang and Alexander N. Gorban, Multivariate Gaussian and Student-t process regression for multi-output prediction, 2019 Springer Science and Business Media LLC

[18]
Duvenaud David, Automatic model construction with Gaussian processes, 2014 PhD Thesis

[19]
Sebastient Petit, Julien Bect, and Emmanuel Vazquez, Parameter selection in Gaussian process interpolation: an empirical study of selection criteria, arXiv Stat.ME

[20]
Mathai Arak and Provost Serge and Haubold Hans, Multivariate Statistical Analysis in the Real and Complex Domains, Chapter 4: The Matrix-Variate Gaussian Distribution, 2022 Springer International Publishing, Pages 217-288

[21]
Marin, Jose Luis and Gaitan, Vicenç and Oms, Guiu and Chiaramello, Marco and Cossart, Quentin and Guironnet, Adrien, An Open Source Tool to Compare Simulators on Large-Scale Cases — Application to Dynawo, 2022 Open Source Modelling and Simulation of Energy Systems (OSMSES)

[22]
Liu, Haitao and Ong, Yew-Soon and Shen, Xiaobo and Cai, Jianfei, When Gaussian Process Meets Big Data: A Review of Scalable GPs, 2020 IEEE Transactions on Neural Networks and Learning Systems

[23]
Virtanen, Pauli and Gommers, Ralf and Oliphant, Travis E. and Haberland, Matt and Reddy, Tyler and Cournapeau, David and Burovski, Evgeni and Peterson, Pearu and Weckesser, Warren and Bright, Jonathan and van der Walt, Stéfan J. and Brett, Matthew and Wilson, Joshua and Millman, K. Jarrod and Mayorov, Nikolay and Nelson, Andrew R. J. and Jones, Eric and Kern, Robert and Larson, Eric and Carey, C J and Polat, İlhan and Feng, Yu and Moore, Eric W. and VanderPlas, Jake and Laxalde, Denis and Perktold, Josef and Cimrman, Robert and Henriksen, Ian and Quintero, E. A. and Harris, Charles R. and Archibald, Anne M. and Ribeiro, Antônio H. and Pedregosa, Fabian and van Mulbregt, Paul and SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, 2020 Nature Methods

[24]
Nowak, Piotr and Hryniewicz, Olgierd, Generalized versions of MV-algebraic central limit theorems, 2015 Kybernetika

[25]
Šošić, Darko and Škokljev, Ivan and Pokimica, Nemanja, Features of Power Transfer Distribution Coefficients in power System Networks, 2014

§ ANNEX A: PROOF OF PROPOSITION 1

The proof relies on the Lyapunov version of the CLT [24]. It states that for $m$ given random independent variables  $X_1,..,X_m$, with expected value $\mu_i$ and variance $\sigma_i^2$, if for some $\delta > 0$, Lyapunov condition is satisfied :

\begin{align*}
s_m^2 &= \sum_{i=1}^m \sigma_i^2 \\
\lim\limits_{\substack{m \to \infty}}& \frac{1}{s_m^2} \sum_{i=1}^m \mathbb{E} \left[ \left| X_i - \mu_i \right|^{2 + \delta} \right] = 0
\end{align*}

\begin{align*}
\frac{1}{s_m} \sum_{i=1}^m (X_i - \mu_i) \longrightarrow_d \mathcal{N}(0, 1)
\end{align*}

In our case, one can easily check that $\forall i \in [1, m]$ :

\begin{align*}
w_i &\sim B\left((2 q_i - 1) \mathbf{p_{safe}} + (1 - q_i)\right) \\
\mathbb{E} \left[ w_i \right] &= (2 q_i - 1) \mathbf{p_{safe}} + (1 - q_i) \\
V(w_i) &= q_i(1-q_i) + \mathbf{p_{safe}} (1 - \mathbf{p_{safe}}) (2 q_i - 1)^2
\end{align*}

Let $A_1,...,A_m$ be $m$ random independant variables following the Bernoulli distributions $B(a_1),...,B(a_m)$, then :

\begin{align*}
s_m^2 &= \sum_{i=1}^m a_i(1-a_i) \\
\mathbb{E} \left[|A_i-\mu_i|^{2+\delta}\right] &= a_i(1-a_i)^{2+\delta} + (1-a_i)a_i^{2+\delta} \\
&\leq 2a_i(1-a_i)
\end{align*}

Since we suppose that the simulator's answer is non-deterministic, there is always a residual uncertainty on the result's prediction. It implies that $q_i(1-q_i)$ does not converge toward $0$ and thus the serie $\sum_{i=1}^m q_i(1-q_i)$ diverges toward $+ \infty$. We deduce that :

\begin{align*}
\lim\limits_{\substack{m \to \infty}} \sum_{i=1}^m \frac{\mathbb{E} \left[|A_i-\mu_i|^{2+\delta}\right]}{s_p^{2+\delta}} &\leq \lim\limits_{\substack{m \to \infty}} \frac{2 \sum_{i=1}^m a_i(1-a_i)}{\left(\sum_{i=1}^m a_i(1-a_i)\right)^{1+\frac{\delta}{2}}} \\
&= 0
\end{align*}

Lyapunov's condition is thus verified for our set of Bernoulli variables. We apply the theorem to $W_1,..., W_m$:

\begin{align*}
&\frac{\sum_{i=1}^m W_i - (2q_i-1)\mathbf{p_{safe}} - (1-q_i)}{\sqrt{\sum_{i=1}^m q_i(1-q_i) + \mathbf{p_{safe}}(1-\mathbf{p_{safe}})(2q_i-1)^2}} \\
= &\frac{\sum_{i=1}^m W_i - (1-q_i) - \mathbf{p_{safe}} \sum_{i=1}^m (2q_i-1)}{\sqrt{\sum_{i=1}^m q_i(1-q_i) + \mathbf{p_{safe}}(1-\mathbf{p_{safe}}) \sum_{i=1}^m (2q_i-1)^2}} \\
= & \sqrt{\frac{\left(\sum_{i=1}^m (2q_i-1)\right)^2}{\sum_{i=1}^m (2q_i-1)^2}} \frac{\frac{\sum_{i=1}^m W_i - (1-q_i)}{\sum_{i=1}^m (2q_i-1)} - \mathbf{p_{safe}}}{\sqrt{\frac{\sum_{i=1}^m q_i(1-q_i)}{\sum_{i=1}^m (2q_i-1)^2} + \mathbf{p_{safe}}(1-\mathbf{p_{safe}})}} \\
= & \sqrt{v_m} \frac{\overline{Z_m}-\mathbf{p_{safe}}}{\sqrt{\sigma_m^2 + \mathbf{p_{safe}}(1-\mathbf{p_{safe}})}}
\end{align*}

§ ANNEX B: PROOF OF PROPOSITION 2

Asymptotically, we know that $\sqrt{v_m} \frac{\overline{Z_m} - \mathbf{p_{safe}}}{\sqrt{\sigma_m^2 + \mathbf{p_{safe}}(1-\mathbf{p_{safe}})}} \sim \mathcal{N}(0, 1)$. To obtain the result, we solve the following inequality :

\begin{align*}
&-q_{\frac{\alpha}{2}} \leq \sqrt{v_m} \frac{\overline{Z_m} - \mathbf{p_{safe}}}{\sqrt{\sigma_m^2 + \mathbf{p_{safe}}(1-\mathbf{p_{safe}})}} \leq q_{\frac{\alpha}{2}} \\
& \implies v_m \frac{\left(\overline{Z_m} - \mathbf{p_{safe}}\right)^2}{\sigma_m^2 + \mathbf{p_{safe}}(1-\mathbf{p_{safe}})} \leq q_{\frac{\alpha}{2}}^2 \\
& \implies \overline{Z_m}^2 + \mathbf{p_{safe}}^2 - 2\mathbf{p_{safe}} \overline{Z_m} \leq \sigma_m^2 \frac{q_{\frac{\alpha}{2}}^2}{v_m} + \mathbf{p_{safe}} \frac{q_{\frac{\alpha}{2}}^2}{v_m} \\
&- \mathbf{p_{safe}}^2 \frac{q_{\frac{\alpha}{2}}^2}{v_m} \\
& \implies \left[ 1 + \frac{q_{\frac{\alpha}{2}}^2}{v_m} \right]\mathbf{p_{safe}}^2 + \left[ -2 \overline{Z_m} - \frac{q_{\frac{\alpha}{2}}^2}{v_m} \right]\mathbf{p_{safe}} \\
&+ \left[ \overline{Z_m}^2 - \sigma_m^2 \frac{q_{\frac{\alpha}{2}}^2}{v_m} \right] \leq 0 \\
\Delta &= \left[ 2 \overline{Z_m} + \frac{q_{\frac{\alpha}{2}}^2}{v_m} \right]^2 - 4 \left[ 1 + \frac{q_{\frac{\alpha}{2}}^2}{v_m} \right] \left[ \overline{Z_m}^2 - \sigma_m^2 \frac{q_{\frac{\alpha}{2}}^2}{v_m} \right] \\
&= 4 \overline{Z_m}^2 + \frac{q_{\frac{\alpha}{2}}^4}{v_m^2} + 4 \overline{Z_m} \frac{q_{\frac{\alpha}{2}}^2}{v_m} - 4 \overline{Z_m}^2 + 4 \sigma_m^2 \frac{q_{\frac{\alpha}{2}}^2}{v_m} \\
&- 4 \overline{Z_m}^2 \frac{q_{\frac{\alpha}{2}}^2}{v_m} + 4 \sigma_m^2  \frac{q_{\frac{\alpha}{2}}^4}{v_m^2} \\
&= 4 \frac{q_{\frac{\alpha}{2}}^2}{v_m} \left[ \overline{Z_m}\left(1-\overline{Z_m}\right) + \sigma_m^2 \left[ 1 + \frac{q_{\frac{\alpha}{2}}^2}{v_m} \right] + \frac{q_{\frac{\alpha}{2}}^2}{4v_m} \right]\\
\end{align*}

We can finally derive the expressions of $p_{min}$ and $p_{max}$ :

\begin{align*}
p_{min} &= \frac{\overline{Z_m} + \frac{q_{\frac{\alpha}{2}}^2}{2v_m} - q_{\frac{\alpha}{2}} \sqrt{ \frac{\overline{Z_m}\left(1-\overline{Z_m}\right)}{v_m} + \frac{\sigma_m^2}{v_m} \left[ 1 + \frac{q_{\frac{\alpha}{2}}^2}{v_m} \right] + \frac{q_{\frac{\alpha}{2}}^2}{4v_m^2}}}{1 + \frac{q_{\frac{\alpha}{2}}^2}{v_m}} \\
p_{max} &= \frac{\overline{Z_m} + \frac{q_{\frac{\alpha}{2}}^2}{2v_m} + q_{\frac{\alpha}{2}} \sqrt{ \frac{\overline{Z_m}\left(1-\overline{Z_m}\right)}{v_m} + \frac{\sigma_m^2}{v_m} \left[ 1 + \frac{q_{\frac{\alpha}{2}}^2}{v_m} \right] + \frac{q_{\frac{\alpha}{2}}^2}{4v_m^2}}}{1 + \frac{q_{\frac{\alpha}{2}}^2}{v_m}} \\
\end{align*}