Optimal Power Flow Based on Physical-Model-Integrated Neural Network
  with Worth-Learning Data Generation
Fast and reliable solvers for optimal power flow (OPF) problems are
attracting surging research interest. As surrogates of physical-model-based OPF
solvers, neural network (NN) solvers can accelerate the solving process.
However, they may be unreliable for ``unseen" inputs when the training dataset
is unrepresentative. Enhancing the representativeness of the training dataset
for NN solvers is indispensable but is not well studied in the literature. To
tackle this challenge, we propose an OPF solver based on a
physical-model-integrated NN with worth-learning data generation. The designed
NN is a combination of a conventional multi-layer perceptron (MLP) and an
OPF-model module, which outputs not only the optimal decision variables of the
OPF problem but also the constraints violation degree. Based on this NN, the
worth-learning data generation method can identify feasible samples that are
not well generalized by the NN. By iteratively applying this method and
including the newly identified worth-learning samples in the training set, the
representativeness of the training set can be significantly enhanced.
Therefore, the solution reliability of the NN solver can be remarkably
improved. Experimental results show that the proposed method leads to an over
50% reduction of constraint violations and optimality loss compared to
conventional NN solvers.
Zuntao Hu, Hongcai Zhang
# Optimal Power Flow Based on

###### Abstract

Fast and reliable solvers for optimal power flow (OPF) problems are attracting surging research interest. As surrogates of physical-model-based OPF solvers, neural network (NN) solvers can accelerate the solving process. However, they may be unreliable for "unseen" inputs when the training dataset is unrepresentative. Enhancing the representativeness of the training dataset for NN solvers is indispensable but is not well studied in the literature. To tackle this challenge, we propose an OPF solver based on a physical-model-integrated NN with worthwhile data generation. The designed NN is a combination of a conventional multi-layer perceptron (MLP) and an OPF-model module, which outputs not only the optimal decision variables of the OPF problem but also the constraints violation degree. Based on this NN, the worth-learning data generation method can identify feasible samples that are not well generalized by the NN. By iteratively applying this method and including the newly identified worth-learning samples in the training set, the representativeness of the training set can be significantly enhanced. Therefore, the solution reliability of the NN solver can be remarkably improved. Experimental results show that the proposed method leads to an over 50% reduction of constraint violations and optimality loss compared to conventional NN solvers.

 Optimal power flow, physical-model-integrated neural network, worth-learning data generation

## I Introduction

Optimal power flow (OPF) is a fundamental but challenging problem for power systems [1]. A typical OPF problem usually involves determining the optimal power dispatch with an objective, e.g., minimizing total generation costs or power loss, while satisfying nonlinear power flow equations and other physical or engineering constraints [2]. Due to the nonlinear interrelation of nodal power injections and voltages, OPF is non-convex, NP-hard, and cannot be solved efficiently [3]. With the increasing integration of renewable generation and flexible demands, uncertainty and volatility have been rising on both the demand and supply sides of modern power systems [4], which requires OPF to be solved more frequently. Thus, fast and reliable OPF solvers have become indispensable to ensure effective operations of modern power systems and have attracted surging interest in academia.

There is a dilemma between the solving efficiency and solution reliability of OPF. Conventionally, OPF is solved by iterative algorithms, such as interior point algorithms, based on explicit physical models [5]. However, these methods may converge to locally optimal solutions. Recently, some researchers have made great progress in designing conic relaxation models for OPF, which are convex and can be efficiently solved [6, 7, 8]. Nevertheless, the exactness of these relaxations may not hold in practical scenarios, and they may obtain infeasible solutions [9]. In addition, the scalability of the conic relaxation of alternating current optimal power flow (AC-OPF) may still be a challenge, particularly in online, combinatorial, and stochastic settings [10].

To overcome the limitation of the aforementioned physical-model-based solvers, some researchers propose surrogate OPF solvers based on neural networks (NNs) [11, 12, 13]. These solvers use NNs to approximate the functional mapping from the operational parameters (e.g., profiles of renewable generation and power demands) to the decision variables (e.g., power dispatch) of OPF. Compared to iterative algorithms, they can introduce significant speedup because an NN is only composed of simple fundamental functions in sequence [12, 13]. However, one of the critical problems of NN solvers is that they may be unreliable if not properly trained, especially for "unseen" inputs in feasible regions due to NNs' mystery generalization mechanism [14].

The generalization of NNs is mainly influenced by their structures, loss functions, and training data. Most published papers propose to enhance the generalization of NN OPF solvers by adjusting the structures and loss functions. Various advanced NN structures rather than conventional fully connected networks are employed to imitate AC-OPF. For example, Overkro _et al._[15] use graph NNs to approximate a given optimal solution. Su _et al._[16] employ a deep belief network to fit the generator's power in OPF. Zhang _et al._[17] construct a convex NN solving DC-OPF to guarantee the generalization of NNs. Jeyaraj _et al._[18] employ a Bayesian regularized deep NN to solve the OPF in DC microgrids. Some researchers design elaborate loss functions that penalize the constraints violation, combine Karush-Kuhn-Tucker conditions, or include derivatives of decision variables to operational parameters. For example, Pan _et al._[11] introduce a penalty term related to the inequality constraints into the loss function. This approach can speed up the computation by up to two orders of magnitude compared to the Gurobi solver, but 18.3% of its solutions are infeasible. Ferdinado _et al._[12] include a Lagrange item in the loss function of NNs. Their method's prediction errors are as low as 0.2%, and its solving speed is faster than DC-OPF by at least two orders of magnitude. Manish _et al._[10] include sensitivity information in the training of NN so that only using about 10% to 25% of training data can attain the same approximation accuracy as methods without sensitivity information. Nellikkath _et al._[19] apply physics-informed
NNs to OPF problems, and their results have higher accuracy than conventional NNs.

The above-mentioned studies have made significant progress in designing elaborate network structures and loss functions. However, little attention has been paid to the training set generation problem. Specifically, they all adopt conventional probability sampling methods to produce datasets for training and testing, such as simple random sampling [10, 11, 12, 13, 15, 17, 20], Monte Carlo simulation [18], or Latin hypercube sampling [19, 16]. These probability sampling methods cannot provide a theoretical guarantee that a generated training set can represent the input space of the OPF problem properly. As a result, probability sampling methods may generate insufficient and unrepresentative training sets, so the trained NN solvers may provide unreliable solutions.

It is important to create a sufficiently representative dataset for training an NN OPF solver. A training set's representativeness depends on its size and distribution in its feasible region [21]. Taking a medium-scale OPF problem as an example, millions of data samples may still be sparse given the high dimension of the NN's inputs (e.g., operational parameters of the OPF problem: renewable generation and power demands at all buses); in addition, because the OPF problem is non-convex, the feasible region of the NN's inputs is a complicated irregular space. Thus, generating a representative training set to cover all the feasible regions of the inputs with an acceptable size is quite challenging. Without a representative training set, it is difficult to guarantee that the NN OPF solver's outputs are reliable, especially given "unseen" inputs in the inference process, as discussed in [22, 23].

To address the above challenge, this study proposes a physical-model-integrated deep NN method with worthwhile data generation to solve AC-OPF problems. To the best of our knowledge, this is the first study that has addressed the representativeness problem of the training dataset for NN OPF solvers. The major contributions of this study are twofold:

1. A novel physical-model-integrated NN is designed for solving the AC-OPF problem. This NN is constructed by a conventional MLP integrating an OPF-model module, which outputs not only the optimal decision variables of the OPF problem but also the violation degree of constraints. By penalizing the latter in the loss function during training, the NN can generate more reliable decision variables.
2. Based on the designed NN, a novel generation method for worth-learning training data is proposed, which can identify samples in the input feasible region that are not well generalized by the previous NN. By iteratively applying this method during the training process, the trained NN gradually generalizes to the whole feasible region. As a result, the generalization and reliability of the proposed NN solver can be significantly enhanced.

Furthermore, comprehensive numerical experiments are conducted, which prove that the proposed method is effective in terms of both reliability and optimality for solving AC-OPF problems with high computational efficiency.

The remainder of this article is organized as follows. Section II provides preliminary models and the motivations behind this study. Section III introduces the proposed method. Section IV details the experiments. Section V concludes this paper.

## II Analysis of Approximating OPF Problems by NN

### _AC-OPF problem_

The AC-OPF problem aims to determine the optimal power dispatch (usually for generators) given specific operating conditions of a power system, e.g., power loads and renewable generation. A typical AC-OPF model can be formulated as

\[\min_{\mathbf{V},\mathbf{S}^{\mathrm{Q}}} \mathrm{C}\,\left(\mathbf{S}^{\mathrm{Q}}\right)\] (1a) s.t.: \[\left[\mathbf{V}\right]\mathbf{Y}_{\text{bus}}^{*}\mathbf{V}^{*}=\mathbf{S}^ {\mathrm{G}}-\mathbf{S}^{\mathrm{L}}, \tag{1b}\] \[\underline{\mathbf{S}}^{\mathrm{G}}\leq\mathbf{S}^{\mathrm{G}}\leq \overline{\mathbf{S}}^{\mathrm{G}},\] (1c) \[\underline{\mathbf{V}}\leq\mathbf{V}\leq\overline{\mathbf{V}},\] (1d) \[\left|\mathbf{Y}_{\text{b}}\mathbf{V}\right|\leq\overline{\mathbf{I}}, \tag{1e}\]

where Eq. (1a) is the objective, e.g., minimizing total generation costs, and Eqs. (1b) to (1e) denote constraints. Symbols \(\mathbf{S}^{\mathrm{G}}\) and \(\mathbf{S}^{\mathrm{L}}\) are \(n\times 1\) vectors representing complex bus injections from generators and loads, respectively, where \(n\) is the number of buses. Symbol \(\mathbf{V}\) is an \(n\times 1\) vector denoting node voltages. Symbol \([.]\) denotes an operator that transforms a vector into a diagonal matrix with the vector elements on the diagonal. Symbol \(\mathbf{Y}_{\text{bus}}\) is a complex \(n\times n\) bus admittance matrix written as \(\mathbf{Y}\) at other sections for convenience. Symbol \(\mathbf{Y}_{\text{b}}\) is a complex \(n_{\text{b}}\times n\) branch admittance matrix, and \(n_{\text{b}}\) is the number of branches. The upper and lower bounds of any variable \(\mathbf{x}\) are represented by \(\bar{\mathbf{x}}\) and \(\underline{\mathbf{x}}\), respectively. Vector \(\overline{\mathbf{I}}\) denotes the current flow limit of branches.

### _AC-OPF mapping from loads to optimal dispatch_

An NN model describes an input-output mapping. Specifically, for an NN model solving the AC-OPF problem shown in Eq. (1), the input is the power demand \(\mathbf{S}^{\mathrm{L}}\), and the output is the optimal generation \(\mathbf{S}^{\mathrm{G}^{\ast}}\). Hence, an NN OPF solver describes the mapping \(\mathbf{S}^{\mathrm{G}^{\ast}}=f^{\text{OPF}}(\mathbf{S}^{\mathrm{L}})\). A well-trained NN should be able to accurately approximate this mapping.

We provide a basic example of a 3-bus system, as shown in Fig. 1, to illustrate how NN works for OPF problems and explain the corresponding challenge for generalization. For simplicity, we assume there is no reactive power in the system and set \(r_{31}=r_{12}=0.01\) ; \(\underline{P}_{i}=0\), \(\overline{P}_{i}=4\), for \(i\in\{1,3\}\), \(P_{2}\in[-7,0]\); \(\underline{V}_{i}=0.95\) and \(\overline{V}_{i}=1.05\), for \(i\in\{1,2,3\}\)

Fig. 1: The 3-bus system.
Then, the OPF model Eq. (1) is reduced to the following quadratic programming:

\[\min_{\mathbf{V},\,\mathbf{P}^{\text{G}}} P_{1}+1.5\times P_{3}\] (2a) s.t.: \[P_{1}=V_{1}\left(V_{1}-V_{2}\right)/0.01+V_{1}\left(V_{1}-V_{3} \right)/0.01, \tag{2b}\] \[P_{2}=V_{2}\left(V_{2}-V_{1}\right)/0.01,\] (2c) \[P_{3}=V_{3}\left(V_{3}-V_{1}\right)/0.01,\] (2d) \[0.95\leq V_{3}\leq 1.05,0.95\leq V_{2}\leq 1.05,\] (2e) \[0\leq P_{1}\leq 4,0\leq P_{3}\leq 4,V_{1}=1, \tag{2f}\]

where \(\mathbf{V}\) is \([V_{1}\,\ V_{2}\,\ V_{3}]^{\top}\), and \(\mathbf{P}^{\text{G}}\) is \([P_{1}\,\ P_{3}]^{\top}\).

Given that \(P_{2}\) ranges from -7 to 0, the 3-bus OPF model can be solved analytically. The closed-form solution of \([P_{1}^{*}\ P_{3}^{*}]\) = \(f_{\text{3-bus}}^{\text{OPF}}(P_{2})\), is formulated as follows:

\[P_{1}^{*}=\left\{\begin{matrix}50-50\sqrt{0.04P_{2}+1},&\text{c1},\\ 4,&\text{c2},\end{matrix}\right. \tag{3}\]

\[P_{3}^{*}=\left\{\begin{matrix}0,&\text{c1},\\ 213\left(1-0.34\sqrt{0.04P_{2}+1}\right)^{2}\\ +50\sqrt{0.04P_{2}+1}-146\end{matrix}\right.,&\text{c2},\end{matrix}\right. \tag{4}\]

where c1 denotes condition 1: \(-3.84\leq P_{2}\leq 0\), and c2 denotes condition 2: \(-7\leq P_{2}<-3.84\).

To further analyze the mapping \(f_{\text{3-bus}}^{\text{OPF}}\), we draw the \([P_{1}^{*}\ P_{3}^{*}]\)-\(P_{2}\) curve according to Eqs. (3) and (4), shown in Fig. 2. Both the \(P_{1}^{*}\)-\(P_{2}\) and \(P_{3}^{*}\)-\(P_{2}\) curves are piecewise nonlinear functions, in which two oblique lines are nonlinear because of the quadratic equality constraints. The reason why the two curves above are piecewise is that the active inequalities change the \([P_{1}^{*}\ P_{3}^{*}]\)-\(P_{2}\) relationship. From an optimization perspective, each active inequality will add a unique equality constraint to the relationship, so the pieces in \(f_{\text{3-bus}}^{\text{OPF}}\) are determined by the sets of active inequalities. In this example, the two pieces in each curve correspond to two sets of active inequalities: \(P_{1}\leq 4\) and \(0\leq P_{3}\). Moreover, the two intersection points are the critical points where these inequalities are just satisfied as equalities.

For a general AC-OPF problem, its input is usually high-dimensional (commonly determined by the number of buses), and its feasible space is partitioned into some distinct regions by different sets of active inequality constraints. From an optimization perspective, a set of active constraints uniquely characterizes the relationship \(\mathbf{S}^{\text{G}^{*}}=f^{\text{OPF}}(\mathbf{S}^{\text{L}})\), and the number of pieces theoretically increases with the number of inequality constraints by exponential order [24, 25, 26]. Therefore, there are massive regions, and each region corresponds to a unique mapping relation, i.e., a piece of mapping function \(f^{\text{OPF}}\).

### _Challenges of fitting OPF mapping by NN_

As shown in Fig. 2(a), to fit the two-dimensional piecewise nonlinear curve of \(f_{\text{3-bus}}^{\text{OPF}}\), we first adopt four data samples by simple random sampling and then use an NN to learn the curve. Obviously, there are significant fitting errors between the fitting and the original lines. Because the training set lacks the samples near the intersections in the curve (where \(p_{2}=-0.384\) in this case), the NN cannot accurately approximate the mapping in the neighboring region of the intersections.

A training set representing the whole input space is a prerequisite for an NN approximating the curve properly. However, it is nontrivial to generate a representative training set by probability sampling. As shown in Fig. 2(a), the intersections of \(f^{\text{OPF}}\) are key points for the representativeness, and the number of intersections increases exponentially with that of the inequality constraints, as analyzed in II-B. When each sample is selected with a small possibility \(\rho\), the generation of a dataset containing all the intersection points are in a low possibility event whose probability is equal to \(\rho^{m}\), where \(m\) is the number of intersections. In practice, the only way to collect sufficient data representing the input space by probability sampling is to expand the dataset as much as possible [27]. This is impractical for large power networks. Therefore, the conventional probability sampling in the literature can hardly produce a representative dataset with a moderate size.

As shown in Fig. 2(b), if we are able to identify the two intersections, i.e., (\(P_{2}=-0.384\), \(P_{1}=4\)) and (\(P_{2}=-0.384\), \(P_{3}=0\)), and include them as new samples in the training dataset, the corresponding large fitting errors of the NN can be eliminated. These samples are termed as the _worth-learning_ data samples. The focus of this study is to propose a _worth-learning data generation method_ that can help identify worth-learning data samples and overcome the aforementioned disadvantage of conventional probability sampling (detailed in the following section).

## III A Physical-Model-Integrated NN with Worth-Learning Data Generation

This section proposes a physical-model-integrated NN with a worth-learning data generation method to solve AC-OPF problems. The proposed NN is a combination of a fully-connected network and a transformed OPF model. It outputs not only the optimal decision variables of the OPF problem but also the violation degree of constraints, which provides guidance for identifying worth-learning data. The worth-learning data generation method creates representative training sets to enhance the generalization of the NN solver.

Fig. 2: Examples of an NN fitting the OPF of the 3-bus system based on (a) simple random sampling, and (b) worth-learning data generation.
### _Framework of the proposed method_

The proposed data generation method has an iterative process, as shown in Fig. 3. First, a training set is initialized by random sampling; second, the physical-model-integrated NN is trained on the training set, where an elaborate loss function is utilized; third, worth-learning data for the current NN are identified; fourth, if worth-learning data are identified, these data are added to the training set and returns to the second step; otherwise, the current NN is output.

The above training process converges until no worth-learning data are identified. This means that the training set is sufficiently representative of the input space of the OPF problem. As a result, the NN trained based on this dataset can generalize to the input feasible set well. The following subsections introduce the proposed method in detail.

### _Physical-model-integrated NN_

In the second step of the proposed method (Fig. 3), the NN is trained to fit the mapping \(\boldsymbol{S}^{\text{G}^{\ast}}=f^{\text{OPF}}(\boldsymbol{S}^{\text{L}})\). To obtain better results, we design a physical-model-integrated NN structure consisting of a _conventional NN module_ and a _physical-model module_, as shown in Fig. 4. The former is a conventional MLP, while the latter is a computational graph transformed from the OPF model.

#### Iii-B1 Conventional NN module

This module first adopts a conventional MLP with learnable parameters to fit the mapping from the \(\boldsymbol{S}^{\text{L}}\) to the optimal decision variable \(\boldsymbol{V}_{\text{NN}}\)[28]. The \(\boldsymbol{V}_{\text{NN}}\) has its box constraint defined in Eq. (1). To ensure that the output \(\boldsymbol{V}_{\text{NN}}\) satisfies this constraint, we design a function \(\text{dRe}()\) to adjust any infeasible output \(\boldsymbol{V}_{\text{NN}}\) into its feasible region, which is formulated as follows:

\[x\leftarrow\text{dRe}(x,\underline{x},\overline{x})=\text{ReLU}(x-\underline {x})-\text{ReLU}(x-\overline{x})+\underline{x}, \tag{5}\]

where \(\text{ReLU}(\text{x})=\text{max}(x,0)\); \(x\) is the input of the function, and its lower and upper bounds are \(\underline{x}\) and \(\overline{x}\), respectively. The diagram of this function is illustrated in Fig. 5.

Applying \(\text{dRe}()\) as the activation function of the last layer of the conventional MLP, the mathematical model of this conventional NN module is formulated as follows:

\[\boldsymbol{V}_{\text{NN}}=\text{MLP}(\boldsymbol{S}^{\text{L}}), \tag{6}\] \[\boldsymbol{V}_{\text{NN}}\leftarrow\text{dRe}(\boldsymbol{V}_{ \text{NN}},\underline{\boldsymbol{V}},\overline{\boldsymbol{V}}), \tag{7}\]

where Eq. (6) describes the conventional model of MLP and Eq. (7) adjusts the output of the MLP.

#### Iii-B2 Physical model module

This module receives \(\boldsymbol{V}_{\text{NN}}\) from the previous module, and then it outputs the optimal power generation \(\boldsymbol{S}^{\text{G}}_{\text{phm}}\) and the corresponding constraints violation \(\boldsymbol{V}\boldsymbol{io}_{\text{phm}}\), where the subscript "phm" denotes the _physical model module_. The first output \(\boldsymbol{S}^{\text{G}}_{\text{phm}}\) is the optimal decision variable of the AC-OPF problem. It can be calculated by \(\boldsymbol{V}_{\text{NN}}\) and \(\boldsymbol{S}^{\text{L}}\), as follows:

\[\boldsymbol{S}^{\text{G}}_{\text{phm}}=[\boldsymbol{V}_{\text{NN}}|\boldsymbol {Y}^{\ast}\boldsymbol{V}_{\text{NN}}^{\ast}+\boldsymbol{S}^{\text{L}}. \tag{8}\]

The second output \(\boldsymbol{V}\boldsymbol{io}_{\text{phm}}\) (termed as violation degree) measures the quality of \(\boldsymbol{S}^{\text{G}}_{\text{phm}}\) and is the key metric to guide the proposed worth-learning data generation (see details in the following subsection III-C). Given \(\boldsymbol{V}_{\text{NN}}\) and \(\boldsymbol{S}^{\text{G}}_{\text{phm}}\), the violations of inequality constraints of the AC-OPF problem \(\boldsymbol{V}\boldsymbol{io}_{\text{phm}}\) are calculated as follows:

\[\boldsymbol{V}\boldsymbol{io}^{\text{S}}_{\text{phm}}=\text{ ReLU}(\boldsymbol{S}^{\text{G}}_{\text{phm}}-\overline{\boldsymbol{S}}^{\text{G}})+ \text{ReLU}(\underline{\boldsymbol{S}}^{\text{G}}-\boldsymbol{S}^{\text{G} }_{\text{phm}}), \tag{9a}\] \[\boldsymbol{V}\boldsymbol{io}^{\text{I}}_{\text{phm}}=\text{ ReLU}(|\boldsymbol{Y}_{\boldsymbol{V}}\boldsymbol{V}_{\text{NN}}|-\overline{\boldsymbol{I}}),\] (9b) \[\boldsymbol{V}\boldsymbol{io}_{\text{phm}}=(\boldsymbol{V} \boldsymbol{io}^{\text{S}}_{\text{phm}}\quad\boldsymbol{V}\boldsymbol{io}^{ \text{I}}_{\text{phm}})^{\top}, \tag{9c}\]

where \(\boldsymbol{V}\boldsymbol{io}^{\text{S}}_{\text{phm}}\) denotes the violation of the upper or lower limit of \(\boldsymbol{S}^{\text{G}}_{\text{phm}}\), and \(\boldsymbol{V}\boldsymbol{io}^{\text{I}}_{\text{phm}}\) represents the violation of branch currents.

_Remark 1_.: _The physical-model-integrated NN is formed by combining the conventional NN module and the physical model module. It inputs \(\boldsymbol{S}^{\text{L}}\) and outputs \(\boldsymbol{S}^{\text{G}}_{\text{phm}}\) and \(\boldsymbol{V}\boldsymbol{io}_{\text{phm}}\), as shown in Fig. 4. Its function is the same as conventional OPF numerical solvers. In addition, it is convenient for users

Fig. 4: The physical-model-integrated NN.

Fig. 5: The dRe() function.

Fig. 3: Framework of the proposed training process.
to directly determine whether the result of the NN OPF solver is acceptable or not based on the violation degree \(\mathbf{Vio_{\text{plm}}}\). In contrast, most NN OPF solvers in the literature are incapable of outputting the violation degree directly [10, 11, 12]._

#### Iii-B3 Loss function

To enhance the training accuracy of the physical-model-integrated NN, we design an elaborate loss function, which consists of \(\mathbf{V}_{\text{NN}}\) from the _conventional NN module_, and \(\mathbf{S^{\text{G}}_{\text{plm}}}\) and \(\mathbf{Vio_{\text{plm}}}\) from the _physical model module_. The formula is as follows:

\[loss=||\hat{\mathbf{V}}-\mathbf{V}_{\text{NN}}||_{1}+||\hat{\mathbf{S}}^{\text{G}}-\mathbf{S }^{\text{G}}_{\text{plm}}||_{1}+\mathbf{Vio_{\text{plm}}}, \tag{10}\]

where \(\hat{\mathbf{V}}\) and \(\hat{\mathbf{S}}^{\text{G}}\) are label values from the training set, which is a ground truth dataset from numerical solvers.

Combining the three terms in the loss function can help enhance fitting precision. As shown in Fig. 6, if the loss function only has the first two items \(||\hat{\mathbf{V}}-\mathbf{V}_{\text{NN}}||_{1}+||\hat{\mathbf{S}}^{\text{G}}-\mathbf{S}^{ \text{G}}_{\text{plm}}||_{1}\) to penalize conventional fitting errors, the predicted value will be in a tiny square space (the red square in Fig. 6) around the label value. From the optimization perspective, the optimal label value is usually on the edge of its feasible region (the blue polyhedron in Fig. 6). This edge through the label value splits the square into two parts: the feasible (blue) part and the infeasible (white) part. Intuitively, we would prefer the predicted values to be in the feasible part. Thus, we also penalize violation degree \(\mathbf{Vio_{\text{plm}}}\) in the loss function to force the predicted values with big \(\mathbf{Vio_{\text{plm}}}\) close to the square's feasible half space for smaller constraint violations.

Although the proposed NN with elaborate loss function has high training accuracy, it is still difficult to guarantee the generalization of the NN OPF solver to the whole input space with conventional random sampling. Therefore, it is indispensable and challenging to obtain a representative training dataset with moderate size to train the proposed NN, which is the focus of the following subsection.

### _Worth-learning data generation_

As shown in Fig. 3, we adopt an iterative process to identify the worth-learning data. For an NN trained in the previous iteration, we utilize its output \(\mathbf{Vio_{\text{plm}}}\) to help identify new data samples that are not yet properly generalized. Specifically, if an input \(\mathbf{S^{\text{L}^{*}}}\) is feasible for the original OPF problem while the current NN outputs a large violation degree \(\mathbf{Vio^{\text{L}}_{\text{plm}}}\), the contradiction means the NN has a large fitting error at \(\mathbf{S^{\text{L}^{*}}}\). This is probably because sample \(\mathbf{S^{\text{L}^{*}}}\) was not included in the previous training set and was not generalized by the NN. Hence, this sample \(\mathbf{S^{\text{L}^{*}}}\) can be regarded as a _worth-learning sample_. Including the sample in the training dataset in the next iteration helps enhance the generalization of the NN.

The key to the proposed worth-learning data generation method is to identify worth-learning samples efficiently. Instead of traversing all of the possible inputs, we maximize \(\mathbf{Vio_{\text{plm}}}\) for a given NN to identify the input with a large violation degree. However, the inputs identified in the maximizing process should be feasible for the original OPF problem. Otherwise, the found inputs might be infeasible and useless for the representation of the training data.

#### Iii-C1 Input feasible set module

To keep the inputs identified in the maximizing process feasible for the original OPF problem, we formulate the _input feasible set module_ to restrict power loads \(\mathbf{S^{\text{L}}}\) to their feasible set. The _feasible set_ is composed of box constraints, current limits, and KCL&KVL constraints, which are transformed from the feasible set of the OPF problem defined in Eq. (1). The partial formulations of the _input feasible set_ are as follows, where the subscript "ifs" denotes the _input feasible set module_:

\[\mathbf{S^{\text{G}}_{\text{ifs}}} =\text{dRe}\left(\mathbf{S^{\prime}}^{\text{G}}_{\text{ifs}},\,\mathbf{ \underline{S}}^{\text{G}},\,\mathbf{\overline{S}}^{\text{G}}\right),\,\,\mathbf{S^{ \prime}}^{\text{G}}_{\text{ifs}}\in\mathbb{R}^{n}, \tag{11a}\] \[\mathbf{V_{\text{ifs}}} =\text{dRe}\left(\mathbf{V^{\prime}}_{\text{ifs}},\,\mathbf{\underline{V} },\,\mathbf{\overline{V}}\right),\,\,\mathbf{V^{\prime}}_{\text{ifs}}\in\mathbb{R}^{n},\] (11b) \[\mathbf{S^{\text{L}}_{\text{ifs}}} =\mathbf{S^{\text{G}}_{\text{ifs}}}\,\,-\,\,\mathbf{[V_{\text{ifs}}]}\mathbf{ Y^{*}}\mathbf{V^{*}}_{\text{ifs}}^{*},\] (11c) \[\mathbf{I_{\text{ifs}}} =\mathbf{Y_{\text{b}}}\mathbf{V}_{\text{ifs}}, \tag{11d}\]

where \(\mathbf{S^{\prime}}^{\text{G}}_{\text{ifs}}\) and \(\mathbf{V^{\prime}}_{\text{ifs}}\) are auxiliary \(n\times 1\) vectors in \(\mathbb{R}^{n}\) and have no physical meaning. Symbols \(\mathbf{S^{\text{G}}_{\text{ifs}}}\) and \(\mathbf{V}_{\text{ifs}}\) are restricted in their box constraints in Eqs. (11a) and (11b). Then the KCL&KVL correlations of \(\mathbf{S^{\text{L}}_{\text{ifs}}}\), \(\mathbf{S^{\text{G}}_{\text{ifs}}}\), and \(\mathbf{V}_{\text{ifs}}\) are described by Eq. (11c). Symbol \(\mathbf{I_{\text{ifs}}}\) in Eq. (11d) denotes the currents at all branches.

The other formulations of the _input feasible set_ aim to calculate \(\mathbf{Vio_{\text{ifs}}}\), the AC-OPF's constraint violations corresponding to \(\mathbf{S^{\text{L}}_{\text{ifs}}}\) and \(\mathbf{I_{\text{ifs}}}\), as follows:

\[\mathbf{Vio^{\text{S}}_{\text{ifs}}} =\text{ReLU}(\mathbf{S^{\text{L}}_{\text{ifs}}}-\mathbf{\overline{S}}^{ \text{L}})+\text{ReLU}(\mathbf{\underline{S}}^{\text{L}}-\mathbf{S^{\text{L}}_{\text{ ifs}}}), \tag{12a}\] \[\mathbf{Vio^{\text{L}}_{\text{ifs}}} =\text{ReLU}(\mathbf{|I_{\text{fis}}|}-\mathbf{\overline{I}}),\] (12b) \[\mathbf{Vio_{\text{ifs}}} =(\mathbf{Vio^{\text{S}}_{\text{ifs}}}\,\,\,\mathbf{Vio^{\text{L}}_{\text{ ifs}}})^{\top}, \tag{12c}\]

where \(\mathbf{Vio^{\text{S}}_{\text{ifs}}}\) denotes the violation of the upper or lower limit of \(\mathbf{S^{\text{L}}_{\text{plm}}}\), and \(\mathbf{Vio^{\text{L}}_{\text{ifs}}}\) denotes the violation of branch current. _Remark 2_.: _This module takes \(\mathbf{S^{\prime}}^{\text{G}}_{\text{ifs}}\) and \(\mathbf{V^{\prime}}_{\text{ifs}}\) as the inputs, and then outputs \(\mathbf{S^{\text{L}}_{\text{ifs}}}\) and \(\mathbf{Vio_{\text{ifs}}}\), as shown in Fig. 7. When \(\mathbf{Vio_{\text{ifs}}}=0\), the corresponding \(\mathbf{S^{\text{L}}_{\text{ifs}}}\) lies in the feasible set of

Fig. 6: Illustration of the effectiveness of the three terms in the loss function.

Fig. 7: The _input feasible set module_.
the AC-OPF problem. To identify feasible \(\mathbf{S}_{\text{ifs}}^{\text{L}}\) in the process of maximizing \(\mathbf{Vio}_{\text{bfm}}\), this module backpropagate the \(\frac{\partial\mathbf{Vio}_{\text{bfm}}}{\partial\mathbf{S}_{\text{bfm}}}\) with \(\mathbf{Vio}_{\text{ifs}}\leq\zeta\) (\(\zeta\) is a small positive tolerance), and then it updates \(\mathbf{S^{\prime}}_{\text{ifs}}^{\text{G}}\) and \(\mathbf{V^{\prime}}_{\text{ifs}}\). As a result, the corresponding \(\mathbf{S^{\prime}}_{\text{ifs}}^{\text{G}}\) is always feasible. Furthermore, because \(\mathbf{S^{\prime}}_{\text{ifs}}^{\text{G}}\) and \(\mathbf{V^{\prime}}_{\text{ifs}}\) are not bounded, changing them can theoretically find any feasible \(\mathbf{S}_{\text{ifs}}^{\text{L}}\)._

#### Iii-C2 Max violation backpropagation

To identify worth-learning data, a novel NN is created by inputting \(\mathbf{S}_{\text{ifs}}^{\text{L}}\) into the physical-model-integrated NN (see Fig. 8). This NN has two outputs, i.e., \(\mathbf{Vio}_{\text{bfm}}\) and \(\mathbf{Vio}_{\text{ifs}}\). The former measures the constraint violation degree of the OPF solution \(\mathbf{S}^{\text{G}^{\text{G}}}\); the latter indicates the feasibility of the OPF input \(\mathbf{S}_{\text{ifs}}^{\text{L}}\). If \(\mathbf{S}_{\text{ifs}}^{\text{L}}\) is a feasible input, i.e., \(\mathbf{Vio}_{\text{ifs}}\leq\zeta\), but the optimal solution \(\mathbf{S}^{\text{G}^{\text{e}}}\) is infeasible, i.e., \(\mathbf{Vio}_{\text{bfm}}\geq\xi\) (\(\xi\) is a threshold), this means the corresponding input is worth learning (i.e., it is not learned or generalized by the current NN). Based on this analysis, we design the loss function \(loss_{\text{max}}\) for max violation backpropagation, as follows:

\[loss_{\text{max}}=\mathbf{Vio}_{\text{bfm}}-\lambda\times\mathbf{Vio}_{\text{ifs}}, \tag{13}\]

where \(\lambda\) is a large, constant weight parameter. When maximizing this loss function, the algorithm tends to find a worth-learning \(\mathbf{S}_{\text{ifs}}^{\text{L}}\) that has small \(\mathbf{Vio}_{\text{ifs}}\) but large \(\mathbf{Vio}_{\text{bfm}}\).

During the max violation backpropagation, the proposed algorithm maximizes \(loss_{\text{max}}\) to update the variables \(\mathbf{S^{\prime}}_{\text{ifs}}^{\text{G}}\) and \(\mathbf{V^{\prime}}_{\text{ifs}}\) by gradient backpropagation until \(loss_{\text{max}}\) converges to the local maximum. After the process, the corresponding \(\mathbf{S}_{\text{ifs}}^{\text{L}}\) is also found. Because the maximizing process can be processed in parallel by the deep learning module PyTorch, the worth-learning samples are found in batch, where the max violation backpropagation uses the previous training set as initial points to identify the new data. Further, the auto-differentiation technique in PyTorch can accelerate the process of parallel computation. Based on these techniques, massive worth-learning data samples are identified efficiently.

### _Overall training process_

The overall training process is presented in Algorithm 1, which first takes an initial training dataset \(\mathbb{D}_{\text{l}}\) (obtained by any conventional sampling method) as input. The learning rate \(\eta\) is equal to \(10^{-3}\), the loss difference tolerance \(\epsilon\) is equal to \(10^{-2}\), the added dataset \(\mathbb{A}\) is empty, and the loss difference \(\Delta L\) is equal to infinity at initialization. The training is performed for a fixed number of epochs (lines 2-5). Then the max violation backpropagation starts to identify worth-learning data (lines 6 and 7) by using the training data as the initial points (line 8) and updating \(\mathbf{S^{\prime}}_{\text{ifs}}^{\text{G}}\) and \(\mathbf{V^{\prime}}_{\text{ifs}}\) until \(\Delta L\) is less than \(\epsilon\) (lines 9-12), which indicates \(loss_{max}\) has converged to the terminal.

```
0:\(\mathbb{D}_{\text{l}}=\left(\mathbf{\hat{S}}^{\text{L}},\hat{\mathbf{V}},\hat{\mathbf{S}} ^{\text{G}}\right)\) Initialization :\(\eta\gets 10^{-3},\epsilon\gets 10^{-2},\mathbb{A}\leftarrow\varnothing, \Delta L\leftarrow\infty\)
1:repeat
2:for epoch \(k=0,1,...\)do
3: Train the NN with \(loss\) Eq. (10):
4:\(w\gets w-\eta\nabla loss\).
5:endfor
6:while\(\Delta L\geq\epsilon\)do
7: Identify data with \(loss_{\text{max}}\) Eq. (13):
8:\(\mathbf{S^{\prime}}_{\text{ifs}}^{\text{G}}\), \(\mathbf{V^{\prime}}_{\text{ifs}}-\mathbf{S}_{\text{ifs}}^{\text{G}}\), \(\mathbf{V_{\text{ifs}}}\leftarrow\mathbf{\hat{S}}^{\text{G}}\), \(\hat{\mathbf{V}}\)
9:\(\mathbf{S^{\prime}}_{\text{ifs}}^{\text{G}}\leftarrow\mathbf{S^{\prime}}_{\text{ifs}}+ \eta\nabla loss_{\text{max}}\)
10:\(\mathbf{V^{\prime}}_{\text{ifs}}=\mathbf{V^{\prime}}_{\text{ifs}}-\eta\nabla loss_{ \text{max}}\)
11:\(\Delta L\leftarrow\left|loss_{\text{max},i}-loss_{\text{max},i-100}\right|\)
12:endwhile
13:\(\{\mathbf{Vio}_{\text{bfm},N}\}\leftarrow f_{\text{filter}}(\mathbf{Vio}_{\text{bfm},N} \geq\xi)\)
14:Collect \(\{\mathbf{S}_{\text{ifs}}^{\text{G}}\}\) corresponding to \(\{\mathbf{Vio}_{\text{bfm},N}\}\) based on the novel NN in Fig. 8
15: Calculate \(\{\hat{\mathbf{V}},\mathbf{S}^{\text{G}}\}\) corresponding to \(\{\mathbf{S}_{\text{ifs}}^{\text{G}}\}\) using numerical solvers
16:\(\mathbb{A}\leftarrow\{\mathbf{S}_{\text{ifs}}^{\text{L}},\hat{\mathbf{V}},\hat{\mathbf{S}} ^{\text{G}}\}\)
17:\(\mathbb{D}_{\text{l}}\leftarrow\mathbb{D}_{\text{l}}\cup\mathbb{A}\)
18:until\(\mathbb{A}\) is \(\varnothing\)
```

**Algorithm 1** Training process of the physical-model-integrated NN OPF solver with worth-learning data generation.

### _Efficiency and convergence of the proposed method_

Unlike general training processes for conventional NNs, the proposed physical-model-integrated NN with worth-learning data generation adopts an iterative training process. It iteratively checks the NN's generalization to the input's feasible space by identifying worth-learning data, as shown in Fig. 3 and Algorithm 1. This difference introduces two critical questions. 1) Efficiency: is the process of identifying worth-learning data computationally efficient? 2) Convergence: is the training set representative of the whole input space after iterations? In terms of the computational efficiency of the proposed method, the theoretical analysis (detailed in the Appendix A) shows it takes no more than 0.08 s to find one sample, which brings little computational burden into the training process. According to the experiment results, the average consumption time for finding one sample is 0.056 s. In terms of the convergence, we prove that the training set would gradually represent the whole input space in the Appendix B, because the number of worth-learning samples identified would converge to zero after a finite number of iterations.

Fig. 8: The novel NN for max violation backpropagation by integrating physical-model-integrated NN with the _input feasible set module_.
## IV Numerical Experiments

The proposed method is evaluated using the IEEE 12-bus, 14-bus, 30-bus, 57-bus, and 118-bus systems. The ground truth datasets are constructed using PANDAPOWER based on a prime-dual interior points algorithm.

### _The efficiency of worth-learning data generation_

As shown in Algorithm 1, the proposed worth-learning data generation (lines 6-12) is the second loop in one iteration (lines 1-18), and the number of initial data points for the generation varies with iterations (lines 8, 15-17). To evaluate the efficiency of the worth-learning data generation, we conduct an experiment on the IEEE 57-bus system in three different iterations to quantitatively measure how much time it takes to finish one worth-learning data generation loop. The time consumption of the data-generation loops in the three different iterations is illustrated in Fig. 9. The x-axis is the number of times the codes are repeated (lines 6-12) divided by 100, which represents the time consumed in one data generation loop; the y-axis is the violation degree. The three lines converge to the terminal stage within 4000 times. The trends are similar: they increase very quickly at first (with 100 epochs) and then approach the local maximum slowly (with 2900-3900 epochs). The inflection points on the three lines are \((1,7228)\), \((1,9065)\), and \((1,5841)\).

In the three iterations, 300, 500, and 800 new data samples are identified. Each data-generation loop in iterations takes 30 s on average to run 3000-4000 times. Hence, one worth-learning data sample costs \((30\times 3)/(300+500+800)\approx 0.056\) s, which introduces little computational burden into the training process compared to the other steps in Algorithm 1. For example, each label value calculated by numerical solvers costs around 1 s (line 14), and the NN training on a dataset with 1100 samples costs around 600 s (lines 2-5). In conclusion, the numerical experiment verifies that the worth-learning data generation brings little computational burden to the training process.

Furthermore, we list the time consumption comparison of the conventional and proposed training processes in Table I, where the conventional training process uses simple random sampling in place of the data generation loop (lines 6-12) in Algorithm 1. By comparing the time consumption of the two methods, we can conclude that the training time of the proposed method only increases by 4%-8%. Hence, these experiments validate that the proposed worth-learning data generation is computationally efficient.

### _Reliability and optimality of the proposed solver_

To validate the superiority of the proposed NN OPF solver (denoted by **Proposed NN**), we compare it with two benchmarks: 1) **B1 NN**, which adopts the conventional loss function and NN model (MLP) with a training dataset generated by simple random sampling; 2) **B2 NN**, which adopts the proposed loss function and physical-model-integrated model with a training dataset generated by simple random sampling.

A particular test set different from the training datasets above is created to examine the effect of these models fairly. The test set has 600 samples that are produced by uniformly sampling 200 points in \([80\%,120\%]\) of the nominal value of one load three times. The other loads in the three times are fixed at light \((80\%\times\text{nominal value})\), nominal \((100\%\times\text{nominal value})\), and heavy \((120\%\times\text{nominal value})\) load conditions. The load sampled has the largest nominal value to cover a big region of the input space. Based on these settings, the test set includes much "unseen" data for those models.

The reliability of the NN OPF solvers is evaluated by the constraint violation degrees on all test data. The optimality loss is evaluated by the relative error between predicted results and label values. For a fair comparison, the three methods all stop their training processes when the value of \(||\hat{\mathbf{V}}-\mathbf{V_{\text{NN}}}||_{1}\) is less than \(2\times 10^{-4}\). In view of the iterative training process, the performance of the three solvers is studied with increasing training data, and the initial NNs are identical because they are trained on an initial dataset with \(N\) samples.

The results are statistically analyzed by creating box plots displayed in Fig. 10. The violation degrees and optimality losses of the results of the NNs from the three methods converge to the terminal stages gradually. The rate of convergence of **Proposed NN** is the largest, that of **B2 NN** is in the middle, and that of **B1 NN** is the smallest.

In Figs. 10(a) to 10(c), the comparison of the last violation degree gives notable results in the three cases. Specifically, the median values in three cases are 7, 15, and 75 for **B1 NN**; 6, 12.5, and 60 for **B2 NN**; and 3.2, 6.1, and 25 for **Proposed NN**, respectively. The novel loss function brings a 19% reduction of violation degree on average by comparing **B1 NN** and **B2 NN**. The proposed training data generation method introduces a 50% reduction of violation degree on average according to the comparison of **B2 NN** and **Proposed NN**. Moreover, the height of the last boxes in each subfigure suggests the robustness of the three solvers, and **Proposed NN**

Fig. 9: Time consumption of the worth-learning data codes in three different iterations. The number of times the sequence codes are repeated in the data generation loop (x-axis) represents the time consumed in one data generation loop; the violation degrees (y-axis) quickly converge to the terminal stage.
has the smallest height in all three cases, which indicates the worth-learning data generation can improve the reliability in encountering "unseen" data from the feasible region.

The comparison of optimality losses is similar to that of violation degrees, as illustrated in Figs. 10(d) to 10(f). The **proposed NN** method has the best results in the three cases, and the final median values of optimality losses are 0.6%, 0.5%, and 0.3% in the three different cases, respectively. The optimality losses of **B2 NN** and **B1 NN** increase by 150%, 66%, and 360% and 142%, 167%, and 460% compared to those of the **proposed NN** method in the three cases.

In conclusion, the proposed physical-model-integrated NN OPF solver with worth-learning data generation can improve the generalization of NN models compared to the conventional NN solvers. Specifically, the proposed method introduces an over 50% reduction of constraint violations and optimality losses in the results on average.

### _Comparison with numerical solvers_

To further evaluate the capability of the proposed method, the next experiment focuses on the comparison with the results of the classical AC-OPF solver based on the prime-dual interior points algorithm and the classical DC-OPF solver with a linear approximation of the power flow equations. The classical AC-OPF solver produces the optimal solutions as the ground truth values, and the DC-OPF solver is a widely used approximation in the power industry. The test set is the same as that in Section IV-B. The performance of the three methods is evaluated by the following metrics: 1) the average consumption time to solve an OPF problem; 2) the average constraint violation degree \(\mathbf{Vio}_{\text{phm}}\), which is calculated by Eqs. (8) and (9) for the two numerical solvers; and 3) the average relative error of dispatch costs. These three metrics are denoted as Time (ms), \(\text{Vio}_{\text{.}}(\text{MW})\), and \(\text{Opt}_{\text{.}}(\%)\), respectively.

The results are tabulated in Table II. The bottom row of the table shows the average results over the three cases. As shown, the proposed method achieves high computational efficiency, which is at least three orders of magnitude faster than the DC-OPF solver and four orders of magnitude faster than the AC-OPF solver. Furthermore, the method also has much lower constraint violations and optimality losses compared with the DC OPF solver. The average Vio. (MW) and Opt. (\(\%\)) of the proposed solver are only 10.882 and 0.462, which are 44% and 18% of those of the DC-OPF solver, respectively.

### _Interpretation of worth-learning data generation_

This subsection interprets why the worth-learning data generated by the proposed method improve the representativeness of the training dataset. The proposed worth-learning data generation method is compared with the conventional simple random sampling method. Without loss of generality, the experiment is conducted on the 14-bus system. Beginning with an identical initial dataset, the conventional and proposed methods generate 100 samples in every step, and there are 8 steps for both. To visualize the representativeness, we draw the distribution of these high-dimensional training samples based on the t-distributed Stochastic Neighbor Embedding algorithm [29, 30], which is a statistical method for visualizing high-dimensional data by giving each data point a location in a two- or three-dimensional map.

The reduced-dimensional data distributions of the conventional and proposed methods are shown in Fig.11. In Fig. 11(a), the data are produced by the simple random sampling method, and their distribution is almost in a "\(\mathbf{\nabla}\)" region, which means the possibility of sampling in this region is high. Furthermore, the new data added in each step overlap with existing data or fill in the intervals. The new data overlapping with existing data are redundant in terms of NN training. The data filling in the intervals may be also redundant when the blanks are generalized well by the trained NN model. In contrast, as shown in Fig. 11(b), the new data generated by

Fig. 10: The violation degree and optimality loss of the results of the NNs trained by three methods change with the number of training data in different cases: (a), (d) IEEE 30-bus; (b), (e) IEEE 57-bus; (c), (f) IEEE 118-bus.
the proposed method in each step hardly overlap with existing data and are usually outside the region covered by the initial data. These new data increase the area covered by the training set so that the training set can have better representativeness of the input feasible region. This explains the effectiveness of the proposed worth-learning data generation method.

## V Conclusion

This study proposes an AC-OPF solver based on a physical-model-integrated NN with worth-learning data generation to produce reliable solutions efficiently. To the best of our knowledge, this is the first study that has addressed the generalization problem of NN OPF solvers regarding the representativeness of training datasets. The physical-model-integrated NN is designed by integrating an MLP and an OPF-model module. This specific structure outputs not only the optimal decision variables of the OPF problem but also the constraint violation degree. Based on this NN, the worth-learning data generation method can identify feasible training samples that are not well generalized by the NN. Accordingly, by iteratively applying this method and including the newly identified worth-learning data samples in the training set, the representativeness of the training set can be significantly enhanced.

The theoretical analysis shows that the method brings little computational burden into the training process and can make the models generalize over the feasible region. Experimental results show that the proposed method leads to over a 50% reduction of both constraint violations and optimality loss compared to conventional NN solvers. Furthermore, the computation speed of the proposed method is three orders of magnitude faster than that of the DC-OPF solver.

## Appendix A Computational Efficiency of Worth-Learning Data Generation

To analyze the computational complexity of the proposed NN model with worth-learning data generation, we adopt a widely used measure--the number of floating-point operations (FLOPs) during the NN model's forward-backward propagation. The total FLOPs of one single layer of a fully-connected NN model can be calculated as follows:

\[\textbf{Forward}:\quad\text{FLOPs}=(2I-1)\times O, \tag{14a}\] \[\textbf{Backward}:\quad\text{FLOPs}=(2I-1)\times O, \tag{14b}\]

where \(I\) is the dimension of the layer's input, and \(O\) is the dimension of its output.

To approximate an OPF mapping based on a 57-bus system, the proposed NN model uses the following structure: \(84\times 1000\times 2560\times 2560\times 5120\times 2000\times 114\). According to Eq. (14), the total FLOPs of the NN per forward-backward process is around \(1\times 10^{8}\). The GPU used in the experiment is the Quadro P6000, and its performance is 12.2 TFLOP/s (\(1\) TFLOP/s \(=10^{12}\) FLOP/s). Using the GPU, we can perform the forward-backward process \(1.22\times 10^{5}\) times per second.

For the worth-learning data generation in Algorithm 1, the forward process is to calculate \(\boldsymbol{V^{i}\alpha_{\text{fs}}}\) and \(\boldsymbol{V^{i}\alpha_{\text{plm}}}\), and the backward process is to update \(\boldsymbol{S^{\prime}\alpha_{\text{fs}}}\) and \(\boldsymbol{V^{\prime}\mathrm{fs}}\) by the gradients. We concatenate \(\boldsymbol{S^{\prime}\alpha_{\text{fs}}}\) and \(\boldsymbol{V^{\prime}\mathrm{fs}}\) as a vector **x**, and we suppose the range of each item in **x** is \([0,10]\), and **x** changes \(10^{-3}\) in each update step. Varying from 0 to 10, it costs \(10^{4}\) times the forward-backward processes. In other words, the algorithm can at least update \(1.22\times 10^{5}/10^{4}\approx 12\) samples in 1 s, so finding one sample costs no longer than 0.08 s.

In practice, there is a slight error between the actual speed in experiments and the theoretical analysis. According to the numerical experiments in Section IV-A, an average of 533 samples are found in 30 s. The average consumption time for identifying one sample is 0.056 s.

Fig. 11: Reduced-dimensional distributions of the training datasets generated by two different methods.
From the analysis presented above, we can conclude that the proposed worth-learning data generation method brings little computational burden into the training process.

## Appendix B Convergence of Worth-Learning Data Generation

This section verifies that the proposed NN with worth-learning data generation can generalize to the whole feasible set. NN models are continuous functions because both linear layers and activation functions are continuous. We define a critical violation value \(\epsilon\) that divides the input space into two types: the covered region (the \(\mathbf{Vio_{\text{phm}}}\) values of all of all of the points are less or equal to \(\epsilon\)) and the uncovered region (the \(\mathbf{Vio_{\text{phm}}}\) values of all of the points are greater than \(\epsilon\)). The boundaries of the two regions consist of the points whose \(\mathbf{Vio_{\text{phm}}}\) values are approximately equal to \(\epsilon\). Using these points as initial points, we can identify points with the local maximum in the uncovered region by max violation backpropagation.

Next, these new points \(\{\mathbf{x}_{1}\}\) (the red points) are added to the training set. After training, the neighborhood of these new points \(\{\mathbf{x}_{1}\}\) would be covered. Due to the generalization of NNs, most points in the area \(\mathbb{S}_{\text{add}}=\{\mathbf{x}|a\times\mathbf{x}_{0}^{\text{ini}}+(1-a)\times \mathbf{x}_{1},0\leq a\leq 1\}\) would also be covered, where \(\{\mathbf{x}_{0}^{\text{ini}}\}\) are the initial points on the boundaries (the black points), as shown in Fig. 12.

Therefore, the area \(\mathbb{S}_{\text{add}}\) is subtracted from the uncovered region. Through iterations, the uncovered region is emptied, and the number of added samples converges to zero.

In practice, we choose the training set instead of the boundary points as initial points for convenience. Although some samples in the training set are not at boundaries, they are eliminated by the filter function, as shown in Algorithm 1. Therefore, the replacement of the boundary points has no impact on the results.
