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

Pierre Houdouin, Manuel Ruiz, Lucas Saludjian, Patrick Panciatici

► To cite this version:

Pierre Houdouin, Manuel Ruiz, Lucas Saludjian, Patrick Panciatici. Certification of MPC-based zonal controller security properties using accuracy-aware machine learning proxies. PSCC, Jun 2024, PARIS, France. hal-04972422

HAL Id: hal-04972422

https://hal.science/hal-04972422v1

Submitted on 1 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# 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, patrick.panciatici}@rte-france.com

***Abstract***—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.

***Index Terms***—Certification of security properties, Congestion management, Multivariate Gaussian processes, NAZA, Proxies

## I. 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 $p_{\text{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 $p_{\text{safe}}$ relies on the law of large numbers. The most basic approach is the brute-force process: scenarios are randomly drawn and simulated until $p_{\text{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 $p_{\text{safe}}$. Using batches of sim-

<figure><img src="image_1.png" /></figure>

<footnote>Submitted to the 23rd Power Systems Computation Conference (PSCC 2024).</footnote>
ulations, 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 p<sub>safe</sub>'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 p<sub>safe</sub> 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.

## II. PRELIMINARIES AND PROBLEM FORMULATION

### A. 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 :

* Z<sub>L</sub> = {L<sub>1</sub>, ..., L<sub>L</sub>} is the set of lines in the considered zone, L is its cardinality

* Z<sub>N</sub> = {N<sub>1</sub>, ..., N<sub>N</sub>} is the set of nodes in the considered zone, N is its cardinality

* P<sub>n</sub> is the allowed renewable power injection at node N<sub>n</sub>

* PA<sub>n</sub> is the available renewable power injection at node N<sub>n</sub>

* P<sub>n</sub><sup>max</sup> is the maximum renewable power that can be produced at node N<sub>n</sub>, and P<sub>n</sub><sup>max</sup> / P<sub>n</sub><sup>max</sup> ∈ [0, 1] is the relative available renewable power injection

* F<sub>l</sub> and F&#773;<sub>l</sub> are respectively the power flow and the IST of line L<sub>l</sub>

### B. Multivariate gaussian process

Consider an unknown function h: ℝ<sup>N</sup> → ℝ<sup>L</sup> that represents the non-deterministic answer y ∈ ℝ<sup>L</sup> of a system to an input x ∈ ℝ<sup>N</sup>. We model y = h(x) + N'(0, σ<sub>0</sub><sup>2</sup>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: μ: ℝ<sup>N</sup> → ℝ<sup>L</sup>

* Its correlation (kernel) function: k: ℝ<sup>N</sup> × ℝ<sup>N</sup> → ℝ

* Its output covariance matrix: Ω, a L × L matrix

Prior knowledge of the system's dynamic can be incorporated into these three elements. For a given x ∈ ℝ<sup>N</sup>, μ(x) represents the prior expected value of h(x). Usually, as no specific information about h is known, the choice is μ = 0, which is also our choice here. The kernel function k(x<sub>1</sub>, x<sub>2</sub>) measures to what extent y<sub>1</sub> and y<sub>2</sub>, the outputs of x<sub>1</sub> and x<sub>2</sub>, are correlated. Its selection has to be adapted to the problem. For our problem, we choose the stationary squared exponential kernel function:

$$k(x_1, x_2) = \sigma_f^2 \exp \left( - \frac{\|x_1 - x_2\|^2}{2l^2} \right) + \sigma_0^2 \delta(x_1, x_2)$$

It is a parametric function with parameters θ = [σ<sub>0</sub><sup>2</sup>, σ<sub>f</sub><sup>2</sup>, l]. Parameter σ<sub>0</sub><sup>2</sup> represents the variance of the system's answer to the same input, σ<sub>f</sub><sup>2</sup> 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 Ω = I<sub>L</sub> implies that coordinates are decorrelated. The MGP then boils down to training L separate univariate GP to predict each coordinate independently. In this work, Ω 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 Ω across the process's iterations is planned. We now define the following variables:

* X<sub>m</sub> = [x<sub>1</sub>, ..., x<sub>m</sub>] the inputs vectors

* Y<sub>m</sub> = [y<sub>1</sub><sup>T</sup>, ..., y<sub>m</sub><sup>T</sup>] the corresponding simulated outputs

* M<sub>m</sub> = [μ(x<sub>1</sub>), ..., μ(x<sub>m</sub>)] the corresponding mean vectors

* Σ<sub>m</sub> = (k(x<sub>i</sub>, x<sub>j</sub>))<sub>i,j∈[1,m]</sub> the covariance matrix between the input vectors

* x<sub>m+1</sub> a new input vector

* Σ<sub>[1,m],m+1</sub> = [k(x<sub>1</sub>, x<sub>m+1</sub>), ..., k(x<sub>m</sub>, x<sub>m+1</sub>)] the covariance vector between the new input and the old ones

The joint distribution of Y<sub>m+1</sub> is a matrix-variate gaussian distribution [20] M N(M<sub>m+1</sub>, Σ<sub>m+1</sub>, Ω) [17]. Similarly to

<figure><img src="image_2.png" /></figure>
the univariate GP, we can derive the conditional distribution of $y_{m+1}|y_1, ..., y_m$:

$$ 
\begin{aligned}
\mu_* &= \mu(x_{m+1}) + \sum_{[1,m],m+1} \Sigma_m^{-1} (Y_m - M_m) \\
\sigma_* &= k(x_{m+1}, x_{m+1}) - \sum_{[1,m],m+1} \Sigma_m^{-1} \sum_{[1,m],m+1} \Sigma_m^{-1} \\
y_{m+1} &\in \mathbb{R}^L \sim \mathcal{N}(\mu_*, \sigma_* \Omega)
\end{aligned}
$$

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

### C. Problem formulation

Given:

* A zone of the power network

* A black-box simulator $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 $p_{\text{safe}}$ to draw a scenario $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:

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

2) Each generator produces the maximum available or allowed power

3) The loads are constants

4) DC approximation: we neglect the reactive power

### Renewable energy generation scenario

We model a renewable energy production scenario $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 $x$ is generated using a truncated multivariate Gaussian distribution.

$$ x = \begin{pmatrix} x_1 \\ \vdots \\ x_N \end{pmatrix} \sim \mathcal{N}_{|\mathcal{X}|} (0, \Sigma_{\text{nodes}}) \rightarrow \begin{pmatrix} PA_1 = x_1 \times P_1^{\text{max}} \\ \vdots \\ PA_N = x_N \times P_N^{\text{max}} \end{pmatrix} $$

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 $S$. It takes as input a renewable power production scenario $x$ and returns the vector of the lines' relative power flow $y \in \mathbb{R}^L$. For line $l$, the relative power flow is $y^{(l)} = \frac{F_l}{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 $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:

$$ f(y) = z = \begin{cases} 1 & \text{if } \forall l \in [1, L], y^{(l)} \le 1 \\ 0 & \text{if } \exists l \in [1, L], y^{(l)} > 1 \end{cases} $$

<figure><img src="image_3.png" /><figcaption>Fig. 1. Black-box function taking as input a renewable energy generation scenario $x$ and returning 0 if a congestion occurred, 1 if security was ensured</figcaption></figure>

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: $p_{\text{safe}} = \mathbb{P}(f \circ S(x) = 1)$. The following section will detail two processes to estimate $p_{\text{safe}}$.

## III. PROCESSES PRESENTATION

### A. Brute-force process

The brute-force process is the vanilla way to estimate $p_{\text{safe}}$. We randomly sample scenarios in $\mathcal{X}$, compute $f \circ S$ using the network simulator, and use the CLT to build a confidence interval around $p_{\text{safe}}$. We iterate the sampling until a satisfying confidence interval is reached. Given the $m$ simulation results

<figure><img src="image_4.png" /></figure>
**Algorithm 1** Brute-force process

**Inputs:**

* scenario sampler $\mathcal{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]$

2. Initialize the bag of outputs $\mathcal{B} = \{\}$

3. **while** not $CI_{good}$(I_m) **do**

4. Draw a scenario $\mathbf{x}_m$ using $\mathcal{X}$

5. Increase $m = m + 1$

6. Compute the outcome $z_m = f \circ \mathcal{S}(\mathbf{x}_m)$

7. Add outcome $z_m$ to the bag $\mathcal{B} = \mathcal{B} \cup \{z_m\}$

8. Compute $I_m$ of level $1 - \alpha$ using $\mathcal{B}$

9. **end while**

$z_1, \dots, 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:

$$I_{\alpha} = [p_{min}, p_{max}]$$

$$p_{min} = \frac{\overline{z_m} + \frac{Q_{\frac{\alpha}{2}}}{2m} - Q_{\frac{\alpha}{2}} \sqrt{\frac{\overline{z_m}(1-\overline{z_m})}{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}(1-\overline{z_m})}{m} + \frac{Q_{\frac{\alpha}{2}}^2}{4m^2}}}{1 + \frac{Q_{\frac{\alpha}{2}}^2}{m}}$$

The brute-force process is summarized in Algorithm 1.

## 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 $\mathbf{y} = \mathcal{S}(\mathbf{x})$ of the simulator to a scenario. We train an MGP-based proxy $\tilde{\mathcal{S}}$ with the previously simulated scenarios.

In this new process, we sample batches of scenarios. For each input scenario $\mathbf{x}_i$, the MGP-based proxy provides the probability density function (PDF) $g_i$ of the output $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 $\mathbf{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

<figure><img src="image_5.png" /><figcaption>Fig. 2. 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</figcaption></figure>

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.

<figure><img src="image_6.png" /><figcaption>Fig. 3. Proxy use decision mechanism</figcaption></figure>

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 \le p_i \le 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.

## C. Central Limit Theorem adaptation

Let $\mathbf{x}_1, \dots, \mathbf{x}_m \in \mathcal{X}$ be $m$ scenarios:

* $z_1, \dots, 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, \dots, w_m \in \{0, 1\}$ denote the answers' prediction

<figure><img src="image_7.png" /></figure>
**Algorithm 2 Proxy-based process**

**Inputs**:

* scenario sampler $\mathcal{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]$

2: Initialize proxy $\tilde{\mathcal{S}}$

3: **while** not $CI_{good}(I_m)$ **do**

4: Draw $m_{batch}$ scenarios $\mathbf{x}_1, ..., \mathbf{x}_{m_{batch}}$ using $\mathcal{X}$

5: Increase $m = m + 1$

6: **for** $i \in [1, m_{batch}]$ **do**

7: Compute the probability $p_i = \mathbb{P}(z_i = 1) = \tilde{\mathcal{S}}(\mathbf{x}_i)$

8: **if** $p_i \le P_{min}$ or $p_i \ge P_{max}$ **then**

9: Draw $w_i \sim B(p_i)$

10: **if** $w_i = 1$ **then**

11: Set $q_i = p_i$

12: **else**

13: Set $q_i = 1 - p_i$

14: **end if**

15: **else**

16: Set $w_i = S(\mathbf{x}_i)$

17: Set $q_i = 1$

18: **end if**

19: Add outcome $(w_i, q_i)$ to the bag $\mathcal{B} = \mathcal{B} \cup \{(w_i, q_i)\}$

20: **end for**

21: Update MGP parameters with simulated samples in $\mathcal{B}$

22: Compute $I_m$ of level $1 - \alpha$ using $\mathcal{B}$

23: **end while**

* 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:

$$ \overline{Z_m} = \frac{\sum_{i=1}^{m} w_i - (1 - q_i)}{\sum_{i=1}^{m} 2q_i - 1} $$

$$ v_m = \frac{(\sum_{i=1}^{m} 2q_i - 1)^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} $$

Then,

$$ \sqrt{v_m} \frac{\overline{Z_m} - \mathbf{p}_{safe}}{\sqrt{\sigma_m^2 + \mathbf{p}_{safe}(1 - \mathbf{p}_{safe})}} \to_d \mathcal{N}(0, 1) $$

**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 | \text{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 $\bar{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:

$$ p_{min} = \frac{\overline{Z_m} + \frac{Q_\alpha^2}{2v_m} - Q_\alpha \sqrt{\frac{\overline{Z_m}(1-\overline{Z_m})}{v_m} + \frac{\sigma_m^2}{v_m} \left[ 1 + \frac{Q_\alpha^2}{v_m} \right] + \frac{Q_\alpha^2}{4v_m^2}}}{1 + \frac{Q_\alpha^2}{v_m}} $$

$$ p_{max} = \frac{\overline{Z_m} + \frac{Q_\alpha^2}{2v_m} + Q_\alpha \sqrt{\frac{\overline{Z_m}(1-\overline{Z_m})}{v_m} + \frac{\sigma_m^2}{v_m} \left[ 1 + \frac{Q_\alpha^2}{v_m} \right] + \frac{Q_\alpha^2}{4v_m^2}}}{1 + \frac{Q_\alpha^2}{v_m}} $$

**Proof**: See Annex B

## IV. COMPUTATIONAL RESULTS

### A. 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 $\mathbf{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

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 4. Selection of the curtailment value $\beta$ for a given scenario based on historical curtailment decisions</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 5. Integral evaluation time depending on the dimension, 10000 averaged simulations per point</figcaption></figure>

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.

$$F_l = \beta \sum_{n=1}^{N} (PTDF_{l,n} + \epsilon_n)x_n$$

$\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 $x_1, ..., x_J$, we compute the minimum required curtailment that avoids congestion: $\beta_j = \arg \max_{\beta \in [0,1]} \{ \beta \sum_{n=1}^{N} (PTDF_{l,n} + \epsilon_n) x_n^j < 1 \}$. Then, for another scenario $x$, we randomly draw its curtailment value $\beta$ among the set of candidates: $Cand(x) = \{ \beta_j, j \in [1, J], |\sum_{n=1}^{N} x_n - x_j^j| < \eta \}$, 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.

### 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 $x_0$.

* At most $N^*$ scenarios

* Only neighbors scenarios that are at a distance less than $d^*: d(x, 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.

### B. 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.

In the beginning, the MGP has no observations from which to learn. The variance of the posterior distribution is

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 6. Evolution of the 200 mobile average of $\sigma_*$ during the proxy-based process</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 7. Evolution of the 200 mobile average of the proxy's prediction entropy during the proxy-based process</figcaption></figure>

$\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 $\Sigma_{[1,m],m+1}^T \Sigma_m^{-1} \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.

In the beginning, the posterior distribution is $N(0, \Omega)$. $\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.

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

<figure><img src="image_14.png" /><figcaption>Fig. 8. Evolution of the total % of simulated scenarios since the beginning during the proxy-based process</figcaption></figure>

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.

### C. 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 $p_{\text{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

<figure><img src="image_15.png" /></figure>
<table><caption>TABLE I
PROCESSES PERFORMANCES COMPARISON</caption><thead><tr><th></th><th>P<sub>safe</sub> = 0.9044</th><th>Proxy-based</th></tr><tr><th></th><th>Brute-force</th><th></th></tr></thead><tbody><tr><th>Number of iterations</th><td>125575</td><td>100000</td></tr><tr><th>p<sub>min</sub></th><td>0.9018</td><td>0.9023</td></tr><tr><th>p<sub>max</sub></th><td>0.9067</td><td>0.9064</td></tr><tr><th>Relative error</th><td>0.05%</td><td>0.041%</td></tr><tr><th>Confidence interval length</th><td>0.0049</td><td>0.0041</td></tr></tbody></table>

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.

<figure><img src="image_16.png" /><figcaption>Fig. 9. Confidence intervals obtained by both processes</figcaption></figure>

## V. 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.

## REFERENCES

[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 (ICAPS)

[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] Sebastien 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] Marín, 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)

<figure><img src="image_17.png" /></figure>
[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, Ilhan 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 :

$$s_m^2 = \sum_{i=1}^{m} \sigma_i^2 \\ \lim_{m \to \infty} s_m^2 \sum_{i=1}^{m} \mathbb{E} [|X_i - \mu_i|^{2+\delta}] = 0$$

Then,

$$\frac{1}{s_m} \sum_{i=1}^{m} (X_i - \mu_i) \to_d \mathcal{N}(0, 1)$$

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

$$\begin{aligned} w_i &\sim B((2q_i - 1)\mathbf{p}_{\text{safe}} + (1 - q_i)) \\ \mathbb{E}[w_i] &= (2q_i - 1)\mathbf{p}_{\text{safe}} + (1 - q_i) \\ V(w_i) &= q_i(1 - q_i) + \mathbf{p}_{\text{safe}}(1 - \mathbf{p}_{\text{safe}})(2q_i - 1)^2 \end{aligned}$$

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

$$\begin{aligned} s_m^2 &= \sum_{i=1}^{m} a_i(1 - a_i) \\ \mathbb{E}[|A_i - \mu_i|^{2+\delta}] &= a_i(1 - a_i)^{2+\delta} + (1 - a_i)a_i^{2+\delta} \\ &\le 2a_i(1 - a_i) \end{aligned}$$

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 series $\sum_{i=1}^{m} q_i(1 - q_i)$ diverges toward $+\infty$. We deduce that :

$$\begin{aligned} \lim_{m \to \infty} \sum_{i=1}^{m} \frac{\mathbb{E}[|A_i - \mu_i|^{2+\delta}]}{s_p^{2+\delta}} &\le \lim_{m \to \infty} \frac{2 \sum_{i=1}^{m} a_i (1 - a_i)}{(\sum_{i=1}^{m} a_i (1 - a_i))^{1+\frac{\delta}{2}}} \\ &= 0 \end{aligned}$$

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

$$\begin{aligned} & \frac{\sum_{i=1}^{m} W_i - (2q_i - 1)\mathbf{p}_{\text{safe}} - (1 - q_i)}{\sqrt{\sum_{i=1}^{m} q_i(1 - q_i) + \mathbf{p}_{\text{safe}}(1 - \mathbf{p}_{\text{safe}})(2q_i - 1)^2}} \\ &= \frac{\sum_{i=1}^{m} W_i - (1 - q_i) - \mathbf{p}_{\text{safe}} \sum_{i=1}^{m} (2q_i - 1)}{\sqrt{\sum_{i=1}^{m} q_i(1 - q_i) + \mathbf{p}_{\text{safe}}(1 - \mathbf{p}_{\text{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{\sum_{i=1}^{m} W_i - (1 - q_i)}{\sum_{i=1}^{m} (2q_i - 1)^2} + \mathbf{p}_{\text{safe}}(1 - \mathbf{p}_{\text{safe}})} \\ &= \sqrt{v_m} \frac{\overline{Z_m} - \mathbf{p}_{\text{safe}}}{\sqrt{\sigma_m^2 + \mathbf{p}_{\text{safe}}(1 - \mathbf{p}_{\text{safe}})}} \end{aligned}$$

ANNEX B: PROOF OF PROPOSITION 2

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

$$\begin{aligned} & -q_{\frac{\alpha}{2}} \le \sqrt{v_m} \frac{\overline{Z_m} - \mathbf{p}_{\text{safe}}}{\sqrt{\sigma_m^2 + \mathbf{p}_{\text{safe}}(1 - \mathbf{p}_{\text{safe}})}} \le q_{\frac{\alpha}{2}} \\ & \implies v_m \frac{(\overline{Z_m} - \mathbf{p}_{\text{safe}})^2}{\sigma_m^2 + \mathbf{p}_{\text{safe}}(1 - \mathbf{p}_{\text{safe}})} \le q_{\frac{\alpha}{2}}^2 \\ & \implies \overline{Z_m}^2 + \mathbf{p}_{\text{safe}}^2 - 2\mathbf{p}_{\text{safe}}\overline{Z_m} \le \sigma_m^2 \frac{q_{\frac{\alpha}{2}}^2}{v_m} + \mathbf{p}_{\text{safe}} \frac{q_{\frac{\alpha}{2}}^2}{v_m} \\ & \qquad - \mathbf{p}_{\text{safe}}^2 \frac{q_{\frac{\alpha}{2}}^2}{v_m} \\ & \implies \left[1 + \frac{q_{\frac{\alpha}{2}}^2}{v_m}\right] \mathbf{p}_{\text{safe}}^2 + \left[-2\overline{Z_m} - \frac{q_{\frac{\alpha}{2}}^2}{v_m}\right] \mathbf{p}_{\text{safe}} \\ & \qquad + \left[\overline{Z_m}^2 - \sigma_m^2 \frac{q_{\frac{\alpha}{2}}^2}{v_m}\right] \le 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} \\ & \qquad - 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} (1 - \overline{Z_m}) + \sigma_m^2 \left[1 + \frac{q_{\frac{\alpha}{2}}^2}{v_m}\right] + \frac{q_{\frac{\alpha}{2}}^2}{4v_m}\right] \end{aligned}$$

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

<figure><img src="image_18.png" /></figure>
$$p_{min} = \frac{\overline{Z_m} + \frac{q_{\frac{\alpha}{2}}^2}{2v_m} - q_{\frac{\alpha}{2}} \sqrt{\frac{\overline{Z_m}(1-\overline{Z_m})}{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}(1-\overline{Z_m})}{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}}$$

<figure><img src="image_19.png" /></figure>