{'fraction_non_alphanumeric': 0.03182665881429132, 'fraction_numerical': 0.006919905771495877, 'mean_word_length': 4.933677926616191, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Closing the Loop: A Framework for Trustworthy Machine Learning in Power Systems

Deep decarbonization of the energy sector will require massive penetration of stochastic renewable energy resources and an enormous amount of grid asset coordination; this represents a challenging paradigm for the power system operators who are tasked with maintaining grid stability and security in the face of such changes. With its ability to learn from complex datasets and provide predictive solutions on fast timescales, machine learning (ML) is well-posed to help overcome these challenges as power systems transform in the coming decades. In this work, we outline five key challenges (dataset generation, data pre-processing, model training, model assessment, and model embedding) associated with building trustworthy ML models which learn from physics-based simulation data. We then demonstrate how linking together individual modules, each of which overcomes a respective challenge, at sequential stages in the machine learning pipeline can help enhance the overall performance of the training process. In particular, we implement methods that connect different elements of the learning pipeline through feedback, thus"closing the loop"between model training, performance assessments, and re-training. We demonstrate the effectiveness of this framework, its constituent modules, and its feedback connections by learning the N-1 small-signal stability margin associated with a detailed model of a proposed North Sea Wind Power Hub system.


I. INTRODUCTION
In order to meet the carbon emission reduction goals set out in the Paris Climate Agreement, governments across the world are targeting massive renewable energy integration projects. Over the next decade, a variety of renewable energy projects related to offshore energy islands [1], GW-scale battery energy storage systems [2], and multi-GW-scale solar parks [3] are set to fundamentally transform global power systems. These large-scale projects will be simultaneously complemented by the deployment of millions of distributed energy resources (DERs). The integration of these renewable energy projects into the existing bulk power grid is ultimately emblematic of the ongoing transition towards zero-carbon energy systems. As legacy grid infrastructure is enriched with these previously unseen yet massively deployed components, system operators are faced with entirely new operational challenges. Ensuring that security constraints, such as N-1 security, smallsignal and voltage stability, and many others, are met in emerging grids, especially when embedding these as constraints in optimal dispatch problems, will quickly become computationally intractable. Examples of increasing intractability, which we find throughout the emerging power grid, necessitate the development of new analysis tools.
With its ability to efficiently learn from complex datasets, provide universal function approximation, and generate predictive solutions on fast timescales, Machine Learning (ML) is emerging as one such tool which can help overcome intractability. Furthermore, emerging research [4], [5] is convincingly showing that embedding surrogate ML models directly inside safety-critical grid routines can transform intractable optimization problems into tractable ones. The recent increase of research on ML within the power systems community is enormous, as comprehensively reviewed by [6]. Several works have already shown ML's potential in real-world applications [7], [8]; however, it is clear that ML has not yet proven ready to be widely adopted in power systems. As [6] elaborates on, there are many existing challenges to be addressed; one of the most critical ones is the building of trust between domain experts and the ML modeling strategies produced by the research community. Trust will be built as the usefulness, efficiency, and reliability of ML models are tested and proven over time.
The barrier associated with building trustworthy models naturally arises from the general complexity associated with ML. Each stage of the ML process, from data collection to model assessment, can be incredibly complex and thus pose a distinct challenge. In this paper, we consider five core challenges associated with building useful and trustworthy ML models. We focus specifically on building ML models learned from physics-based simulations in the context of safety-critical applications, e.g., security assessment and operational decision making. After stating the challenges, we then identify suitable solutions from the literature to solve these problems. The individual solutions are then organized into a coherent framework. Fig. 2 shows an overview of this framework, illustrating the various steps (subsequently called modules) that are involved in building and utilizing a ML model. For each module, substantial bodies of research already exist; some of these works are tailored for power system applications, but most are not application specific.
Combining all of these modules into a single ML pipeline forms a process that can be difficult to handle. However, in order to produce ML models which are useful, efficient, and reliable, we believe that the following two approaches are both vital: 1) the building up of a coherent framework which considers all aspects of ML, and 2) the designing of feedback loops within the framework in order to create synergies between the modules. Fig. 2 illustrates these two approaches by 1) considering the entire ML pipeline and 2) adding additional loops in the process signified by the red arrows to enhance model performance. The overall goal is to make each module more robust and efficient, thus enhancing the overall process itself.
The core contributions of this paper follow: • We consider five core challenges associated with using ML to learn from physics-based simulation data, and we present a coherent modular framework which utilizes state-of-the-art solutions methodologies to tackle the identified challenges. • By utilizing feedback mechanisms, we identify how "closing the loop" can help indicate areas of weakness within a ML model and how targeted data re-sampling and model re-training can greatly benefit the model. • We explicitly demonstrate how the ML framework can be used to learn security-constrained stability margins from a highly detailed dynamical simulation model (i.e., a proposed North Sea Wind Power Hub system). • Finally, we provide a public, open-source repository on Github [9] so that other researchers within the power systems community can test their domain-specific ML challenges using the proposed fully-integrated framework. The remainder of this paper is structured as follows. In Section II, we provide an overview of five central motivating challenges in ML. In Section III, we provide a literature review to identify state-of-the-art solutions for each of these challenges. We subsequently combine leading methods from the literature into a single coherent framework, and we explain how loop-closing can help link these methods together. In Section IV, we present the technical details associated with the chosen case study: learning a security constrained smallsignal stability margin associated with a North Sea Wind Power Hub model. We then demonstrate and compare the improved performance of the resulting models in Section V. In Section VI, we elaborate on potential routes forward for uses of ML within the power systems research community, and we offer conclusions in Section VII.

II. MACHINE LEARNING CHALLENGES
The following subsections provide an overview of the most pressing challenges associated with building trustworthy ML models. In particular, we focus on the challenges which naturally arise when deploying a "fully integrated" ML framework, i.e., a framework that spans the full spectrum of a ML model's intended training and usage, as displayed in Fig. 2. The five key topics we focus on hereafter are summarized below: • Dataset creation: collection of data for the modelling process, e.g., through measurements or simulation; • Data pre-processing: feature engineering of the dataset; • Model training: fitting a parameterized ML model to the provided dataset in a supervised fashion;

A. Dataset creation
The starting point for every supervised ML algorithm is a dataset, either collected from measurements or generated from simulations. In this paper, we focus on simulation-based data collection, which presents the opportunity for creating arbitrarily large datasets. Practical limitations arise, however, from the computational burden of creating large and comprehensive datasets. The often cited "curse of dimensionality" quickly renders exhaustive sampling approaches intractable 1 ; this becomes particularly pressing when the underlying physical function exhibits a high degree of non-linear behavior, or if there is a very large number of input features. As the feature number grows, the model requires a combinatorially increasing number of samples to maintain a dense coverage across the input domain. To overcome this, one may wish to avoid "trialand-error" sampling by directly choosing inputs which map to desirable outputs (i.e., ones close to the region of interest). Inverting nonlinear functions (e.g., time domain simulations) from output to input, however, can be highly intractable in itself. A variety of methods have therefore been developed for more efficient sampling. With such methods, particular effort has been devoted to (i) excluding non-informative regions of the input domain on one hand, and (ii) targeting highlyinformative regions more intensely on the other. Finally, for many power system applications, the initial sampling space is often unbalanced [10], with a large volume of input points mapping to outputs which are infeasible or far from the region of training interest [11] (e.g., when learning a stability boundary).

B. Dataset pre-processing
Before starting to learn a ML model, a pre-processing shall ensure that the subsequent learning task is well posed. This primarily relates to understanding and potentially exploiting or countering statistical properties of the dataset. The first step entails a data exploration in which statistical properties such as class imbalance or differences of scales should be identified. By subsequently applying, for example, transformations or feature selection and engineering, the dataset can be brought in a more suitable form. Another conclusion can be the need for an enrichment of the dataset. The final step of the preprocessing is to split the dataset into subsets for training, validation, and testing for the subsequent process. It is important that previously identified statistical characteristics are captured well in each of the three subsets, as the learning and assessment tasks heavily rely on their distributional similarity in order to avoid biases in the process. Only then can the ML models achieve good and reliable performance when deployed in real world application.

C. Model training
To describe a dataset through a model, one must find some function g(·) that maps the input features x to predicted outputsŷ x . The function itself can take a variety of different forms, but its predictions depend on a set of parameters θ. The objective is to minimize the difference between the ground truth y x and prediction by adjusting the parameters θ: This formulation is equally valid for a simple linear regression model and for a highly complex deep neural network (NN); however, such models differ significantly in the degree of complexity they can represent. Finding an adequate level of representation capacity and a suitable set of parameters are the core problems associated with training a NN. An overly simplistic model will capture general trends, but it will miss more detailed aspects of the underlying mapping ("underfitting"). An overly complex model, however, might perform well on training data but poorly on unseen testing data ("overfitting"). To balance this trade-off and control the influence of other parameters such as optimizer settings, or regularization strengths, so-called "hyper-parameter tuning" is conducted. It serves as a model selection process by judging the model performance on the validation dataset. This model selection can require many repeated training runs to identify well-performing setups, and it is essential for ensuring repeatability and robustness. Achieving these two key properties can enable fair comparisons between different ML models.

D. Performance assessment
Data handling practices strongly influence the performance assessment of ML models. If the statistics of the test dataset are known to match the ones of the process we intend to model, we can formulate a valid expectation on the model performance. If the underlying statistics do not match (i.e., the model is deployed "out of sample"), then the trained model can achieve an excellent performance on the test dataset but perform very poorly when deployed. This leads to a mistaken overtrusting of the model. In either case, this assessment is a purely statistical quantity, and it offers no indication of the model's expected behavior at unseen points. This can lead to surprisingly bad predictions and unexpected behavior. Accordingly, such statistical methods cannot offer rigorous guarantees for NN performance across the full input domain. In safety-critical applications, stronger guarantees, such as Neural Network verification, are often desired, potentially across the global operating region. This is one of the topics we will discuss in this paper.

E. Model embedding
Typically, the goal of deploying a trained ML model serves either (i) solving a previously unsolvable or intractable problem, or (ii) estimating a solution much faster than established methods can. The latter can be especially helpful for screening scenarios where evaluation speed is critical. One challenge, however, is that ML models do not necessarily offer insights ACCEPTED FOR PRESENTATION IN 11TH BULK POWER SYSTEMS DYNAMICS AND CONTROL SYMPOSIUM (IREP 2022), JULY 25-30, 2022, BANFF, CANADA 4 as to why a certain solution was predicted. Another general challenge is associated with "embedding" the ML model inside of a broader optimization or control problem. Embedding allows ML models to be flexibly utilized in a diversity of contexts. Such embedding, however, required finding suitable representations to incorporate the NN's predictions into such problems, since NNs in their common form are not well-suited to be solved via standard numerical solvers (e.g., Newtonbased methods). A final challenge relates to ensuring the continued applicability of the learned model in its intended real world usage and taking actions in case problem characteristics changed since the ML model was trained and assessed.

III. OVERCOMING MACHINE LEARNING CHALLENGES: INTEGRATING SOLUTIONS INTO A WORKFLOW
In the following subsections, we review state-of-the-art solutions from the literature associated with the main challenges of the five modules outlined in Section II. We subsequently organize selected solutions into a modular workflow, and we identify the feedback loops which can create productive synergies between the modules. First, though, we introduce useful mathematical notation to properly describe the modeling challenges.
Mapping problem: We define a physical power system model M whose mapping of interest f : with input domain X : x ∈ R n and output range Y : y x ∈ R m . Usually, we are interested in learning the mapping only on some bounded input domain: x ≤ x ≤ x. We assume there is some subset of the output space that has a disproportionately "high information content" [12]; this could be near a stability boundary. We denote this output region as Y * ⊆ Y. Neural network: We generally refer to the NN mapping bŷ The function NN(·) could take many forms, but we assume in the following a simple feed-forward network, with ReLU (Rectifier Linear Unit) activation functions, formulated as where K represents the number of layers, each of which is defined by a weight matrix W k and a bias vector b k . Their sizes are [N k+1 × N k ] and [N k+1 × 1] respectively, where N k is the number of neurons in the k-th layer. The NN training problem then reads as Mixed Integer Linear Program (MILP) reformulation: When using ReLU activation functions (as shown by (6) in our case), it is possible to reformulate (4)-(7) as a MILP. It hinges around replacing the ReLU activation function by a set of linear constraints and N k binary variables [13]: This exact reformulation of the NN prediction to a MILP allows us to integrate it to any optimization or control problem, which we can exploit, e.g., for the model assessment and embedding modules.

A. Dataset creation and enrichment
A number of classical ML textbooks (e.g., [14]) provide a review of data collection practices which are rigorously rooted in statistical sampling theorems. In this module, we assume a dataset can be created by querying a physics-based model M.
That is, we can arbitrarily sample in X ⊆ R n to produce a data point y x from (2). In III-A1, we first present an overview of methods from the literature used for producing representative datasets from physics-based power system models. In III-A2-III-A4, we then present ideas concerning how to "enrich" the dataset, either based on directed walks (DWs), or by using the NN itself to recommend targeted regions of retraining.
1) Sampling approaches in the power systems literature: Using tools from mathematical finance, [15] employs importance sampling to generate probabilistic security assessment risk data. Importance sampling is also employed in [16] for decision tree training, but a two stage procedure is employed in order to bias the sampling routine towards the high information content input space.
To sample quickly from the feasible power flow space, [17] proposed a hypersphere "grid pruning" algorithm rooted in SOCP and Lasserre hierarchies. This work was extended in [11], where SDP-based hyperspheres are constructed in order to classify regions of infeasibility, thus enabling rapid rejection sampling of future input points; this work found power flow solutions that were both feasible and N-1 secure. To improve upon these methods, authors in [10] utilize the quadratic convex (QC) relaxation [18] of the power flow problem in order to generate infeasibility certificates based on separating hyperplanes (rather than hypershperes). Hit and run sampling [19], a Monte Carlo sampling approach for generating asymptotically distributed datasets, has also found application in the power systems literature. Most recently, the developers of the OPF-Learn tool [20] use hit and run sampling to "uniformly and quickly" sample high dimensional polytopes. Using SOCP relaxations, separating hyperplanes are then constructed for infeasibility classification.
Other works related to dynamic and static security assessment have focused more sharply on collecting targeted datasets close to a pre-defined stability boundary (e.g., a smallsignal [11], transient [21], or steady state [22]  addressing dataset imbalances by including the associated risk and cost of a data point in the model learning and assessment [23]. In [12], a dynamic security assessment dataset is iteratively enriched through the addition of points which are at the midpoint between secure and insecure operating points (OPs). To improve upon this method, [11] utilized a so-called "directed walk" method which iteratively pushed operating points close to a small-signal stability boundary via gradient descent.
2) Enhancing datasets with Directed Walks (DWs): When pre-processing datasets for learning, it can become clear that the statistical properties of the dataset are not satisfactory, e.g., due to an insufficient representation of critical areas in the output domain. Therefore, one of the approaches we utilize in this paper is the directed walk methodology developed in [11] in order to target data points which lie close to a pre-defined stability threshold.
3) Enhancing datasets with neural network informed (NI) points: The NN itself, even in a not-fully trained state, can provide an indication of which points in the input domain map to regions in the output domain that are most critical for obtaining a balanced and informative dataset. These mappings might not be perfectly accurate, but they can be used to perform a pre-selection of points that are first re-evaluated using the physical model and then added to the dataset. Algorithm 1, which is first proposed in this paper, describes the essential steps of this approach, which can be viewed as a practical way to close the loop between model training, NN prediction, and model retraining.

Algorithm 1 Network Informed (NI) dataset enrichment
Require: NN, target region Rc, Nsamples 1: Create a large number (e.g., 1'000'000) of random samples x in the input domain X 2: Evaluate neural networkŷx = NN(x) 3: Select all samples x withŷx ∈ Rc − → xc 4: Select Nsamples NI samples out of xc − → xNI 5: Evaluate ground truth yx,NI = f (xNI ) 6: return NI dataset (xNI , yx,NI ) 4) Enhancing datasets using verification informed (VI) points: This approach is also proposed for the first time in this paper, and it also aims to "close the loop" (see Fig. 2). Compared to the NI approach, instead of simply choosing sampled data based on NN evaluations, we now evaluate entire regions of the NN input domain by searching for certain predictions. By utilizing the NN's exact MILP reformulation from (10a)-(10e), it is possible to pose an optimization problem that allows us to test entire regions of the NN for certain behavior (see, e.g., [24], where the MILP reformulation is used to extract worst-case performance guarantees). With this reformulation, we can use optimization-based tools to efficiently determine regions in the input domain that either are, or are not, of interest, and we can find so-called adversarial points (i.e., points whose NN predictions deviate sharply from other points close by). By feeding this information back into the dataset creation module, we can enrich the dataset with more targeted re-sampled data points. Algorithm 2 illustrates the procedure, where we search for values of which certify particular regions of NN prediction (see, e.g., (19) and the associated illustration in Fig. 5).

Algorithm 2 Verification Informed (VI) dataset enrichment
Require: NN, target region Rc, Nsamples 1: Create a large number (e.g., 1'000'000) of random samples x in the input domain X 2: for each reference point x0 do The focus of pre-processing the dataset lies in preparing the dataset in such a way that the subsequent training problem can be solved efficiently. This step can range from applying straightforward best practices to highly involved techniques; in any case, having a profound understanding of the present dataset can prove vital for all subsequent steps. Therefore, it is worth combining pre-processing with data exploration techniques so one can be aware of potential statistical particularities of a dataset.
In terms of practical aspects of pre-processing, the training can usually be aided by standardizing the input data; ideally, the input features are also uncorrelated [25]. The reasons behind this relate to the optimization algorithms used in the training process performing more efficiently. Furthermore, standardization helps to avoid unintended effects due to scaling when using regularization in the training. A more fundamental and theoretical question relates to the selection or engineering of input and output features, i.e., the question of which features or combinations of features in the dataset are the most important ones for accurate predictions. The idea of "engineering" the learning task links to the assumption that the mapping which we attempt to learn could be expressed in a lower dimension if an appropriate set of coordinates, or more formally, the appropriate manifold, is identified. Selecting and engineering input and output features plays an important role for subsequent steps of the process, and it can greatly improve of hurt the chances of success. Refs.
[26]- [28] describe relevant techniques in a general ML context, for reviews for power systems we refer to [6], [29]. In a setup where a dataset is generated from simulations, using the available domain knowledge can render an advanced feature selection process mostly unnecessary; nonetheless, a different representation of the data is worth considering to ease the learning task. Pre-processing can also include the "cleaning" of a dataset. This could entail, among other things, removing or imputing missing values or properly treating of outliers. This delicate procedure is beyond the scope of this work, since we assume a simulation-based dataset creation. However, for measurement-based datasets this "cleaning" is critical.
At the stage of the pre-processing it is also advisable to decide on the methodology of model training, model selection and model assessment. If the base dataset is sufficiently large in terms of data points, the preferred set-up is to perform a random three-fold split of the dataset into training, validation, and testing subsets. Common practice is to retain 20% of the dataset as the hold-out or test dataset to obtain an unbiased estimation on the generalization error. Quoting Hastie et al. [26] for a best practice: "Ideally, the test set should be kept in a 'vault', and be brought out only at the end of the data analysis." Data analysis, in this context, corresponds to the model assessment, i.e., the model that seems to describe the analyzed data "best". The remaining dataset is then further split, for example, at a ratio 80% to 20% for the training and validation set. These numbers should not be taken too literally, but rather act as guides, since they can be problem and setup specific. The respective sections in [14] and in particular [26] provide a good introduction into the underlying concepts and statistical considerations for selecting splitting ratios. The critical aspect in splitting is that the three subsets are distributionally as similar as possible. Otherwise, the model selection and generalization error estimation can be severely biased. For setups with fewer data points, [26] also elaborates on resampling methods such as cross-validation and bootstrapping techniques. Instead of random splitting, advanced tools, such as Algorithmic Splitting [30], are also being developed to determine superior dataset splitting for particular training problems. For setups with access to the data generating process, i.e., the underlying physical model for creating data points, usually the three-fold split is possible.

C. Model training and physics-based regularization
In its vanilla form, the training process is guided purely by the objective of minimizing the loss function L given by where y x stems from a ground truth dataset, andŷ x represents the NN's prediction. Even though it often possible to drive the loss value on the training dataset close to zero, the generalization error, i.e., a measure of the performance on the testing dataset, can turn out to be very high. This describes the phenomenon of "over-fitting" which we try to avoid by controlling the generalization error on the validation set during the training. This can take the form of "early stopping" as explained in [25], where the training is stopped as soon as the the error on the validation set does not improve significantly anymore. Another method on the level of the training is to introduce regularization, such as applying 1norm or 2-norm regularization to the weights of the NN (i.e., penalizing the magnitudes of the weights). Regularization influences the bias-variance-trade-off of the models [26], and if its strength is well-tuned, it helps the training algorithm generate models which are both more generalizable and more parsimonious [31]. This leads directly to a third area of managing generalization error: "hyper-parameter tuning". This refers to the problem of finding suitable parameters for the training process, e.g., optimizer settings or regularization strength, and for the parameters governing the model complexity, e.g., the size of the NN. Since these hyper-parameters influence the model performance in terms of generalization capacity in non-trivial ways, one usually varies them either by random trials, a structured grid search, or by relying on a Bayesian approach. The validation dataset serves once more as a check for the generalization capacity and is used for model selection, where the term "model" also incorporates the training hyperparameters.
As previously mentioned, regularization is a powerful tool for increasing generalization capacity. Since power system problems often fall into domains governed by physical models, the addition of physics-based regularization might be possible and can potentially offer performance improvements. To add this regularization, the loss function L is, for example, extended by a comparison between the derivatives of the physical mapping and NN mapping, as in (12). More generally, regularization can be provided through a system-specific function h(·) which links the inputs, outputs, and potentially their derivatives. Thereby, differential and algebraic relations among the variables of the governing model can be incorporated and deviating predictions get penalized: These additional terms can be weighted by α and β, and they can be used flexibly depending on physical model availability. The NN's derivative can usually be computed using so-called automatic differentiation [32]; for the physical model, gradients can either be derived directly from the model formulation, e.g., for differential equations, or they can be evaluated numerically. 1) Regularization approaches in power systems: Many recent results have focused on a particular regularization architecture, known as Physics Informed Neural Networks (PINNs) [33], where partial/ordinary differential equation (P/ODE) functions explicitly regularize NN training. While this technique is relatively new, it has quickly attracted thousands of followup extensions and associated publications [34].
PINNs were first applied to power systems in [35] to predict swing equation dynamics. They have since additionally been used for system identification [36], transient stability predictions [37], and for learning grid dynamics without simulation data [38]. Beyond ODE simulation and trajectory prediction, physics and sensitivity informed methods have also been utilized for regularizing models related to power distribution grid optimization [39], ACOPF [40], [41], DCOPF [42], parameter estimation [43], and risk-aware voltage optimization via "riskregularization" [44]. In this paper, we utilize numerically computed gradient information in order to directly regularize the training process. Incidentally, much of this gradient information is "left over" from the directed walk dataset enrichment step, so it efficiently serves as a double use.

D. Neural network verification and performance assessment
Verifying the performance of NN-based ML models has become an increasingly important subfield of AI research; an excellent overview is provided in [45]. In this subfield, researchers have predominantly focused in the past on quantifying the so-called adversarial robustness of a NN [45], other approaches attempt to quantify uncertainty through algorithms using conformal predictions as introduced by Vovk et al. [46] and extended in [47]; such methods, however, cannot offer rigorous guarantees for NN performance. In safety-critical applications, stronger guarantees are often desired, potentially across the global operating region.
To generate provable robustness and worst-case performance guarantees, researchers have developed methods which transform the NN mapping into an equivalent set of optimization-based constraints. Such methods have utilized convex relaxations of activation functions and Lagrangians, norm-bounding, and abstract interpretations [48]- [58]. In the field of Reinforcement Learning (RL), researchers have also focused on training agents which are safe by design (i.e., "safe RL") [59], via e.g., learning control policies based on Lyapunov design principles [60] and reachability analysis [61]. Such methods can suffer from (i) learning suboptimal control policies and (ii) an inability to provide guarantees against all possible contingencies [45].
An alternative verification approach relies on modelling the ReLU-based NN as a binary MILP [13], [62]- [64]. Unrelaxed MILP verification is known as a "complete verification algorithm" in the literature, because it utilizes an exhaustive search of the constrained input space and therefore can produce exact, worst-case guarantees [50]. While MILPs are generally NPhard to solve, MILP solvers have become tremendously powerful over the last decade. With light sparsity regularization and pruning, the authors in [64] could solve verification problems on NNs with 1000+ ReLUs in 10s of seconds.
Power system operation and control problems are considered safety-critical applications; therefore, ML models should meet strong verification standards if they are ever to be realistically deployed within these realms. Complete verification methods (reviewed above) can generate the strongest and tightest possible guarantees across large operational possibilities and are, therefore, considered a promising research direction. This approach was first applied to power systems in [65], where the authors computed formal guarantees for the performance of classification NNs in the context of security assessment. This method was extended in [24] and [41], where the authors computed worst-case performance guarantees for a regression NN which was trained to solve the DC-and AC-OPF problems, respectively. A feedback procedure is developed in [42], which finds the worst-case performing input point, and then adds this points (and its true ground-truth output) into the training set.
In this paper, we assume access to the underlying physical model (i.e., the dynamic power system model describing the North Sea Wind Power Hub) for training purposes; however, the mapping from input (power dispatch and control parameters) to output (worst-case security constrained damping ratio) is overly intractable to be used for generating global worst-case performance guarantees. Accordingly, we explore alternative MILP-based verification approaches which iteratively (i) locate regions where training data is lacking and then (ii) solve a bounded worse-case performance problem in order to add targeted points back into the training dataset.

E. Embedding neural networks in optimization problems
Increasingly, ML models have been deployed as surrogate models which replace intractable constraints inside of challenging optimization problems. Generally, this procedure is referred to as "embedding" a ML model inside of an optimization problem [66], and it has enjoyed a growing spectrum of uses, e.g., via MILP reformulation of multivariate adaptive regression splines [67], co-prediction and optimization of unknown parameter values [68], and optimization over tree ensembles [69]. In particular, reformulating ReLU-based NNs into equivalent MILP constraints [13] is becoming a popular tool for embedding NNs inside of optimization problems [70]- [72]. While this approach is highly flexible, the computational expense of solving the associated MILP grows quickly with NN size, becoming nearly intractable for NNs with thousands of nodes [66].
NN embedding has only been applied within the power systems research community in few applications so far. Authors in [73] approximate the nonlinear function between microgrid operating point and frequency nadir in order to encode frequency constraints in a scheduling problem. Surrogate frequency constraints are then encoded into the linearized Unit Commitment problem in [74]. In order to solve an AC Unit Commitment problem, [4] replaces the AC power flow equations with a piecewise linear approximation learned by an optimally compact NN. In [5], a NN learns to classify regions of small-signal stability in order to solve a security constrained OPF problem. Input convex NNs, as popularized in [75], are used in [76] to learn a reactive power control law for optimal voltage regulation. [77] learns grid-following converter dynamics with a NN, and then uses a MILP reformulation to estimate associated critical clearing times.
This paper introduces the benefits that trustworthy ML can offer in power systems, and it demonstrates them on offshore wind hubs in order to support planners and operators who must safely and optimally dispatch their systems. Accordingly, after building a ReLU-based NN which maps from dispatch decisions and control parameters to damping ratios, we embed this mapping inside of an illustrative security-constrained dispatch problem. This embedding process utilizes the MILP reformulation proposed in [13], and it highlights the capacity for a NN to learn intractably complex constraint mappings. We now summarize the core modules in the proposed framework, as depicted by the workflow in Fig. 2.
Usage criteria: This framework is most applicable for building surrogate ML models from physics-based simulation data whose underlying input/output mapping can be arbitrarily (although, potentially expensively) queried. Furthermore, we assume that naively generated training datasets will be unbalanced (i.e., they will not adequately represent the regions of highest interest), and we therefore need to develop approaches to address that. Finally, we assume the resulting ML model will be embedded inside of an optimization or control problem.
Framework summary: We begin by defining bounded input X and output Y domains of interest associated with a physical mapping (2). In particular, the set of output features are chosen to be sufficiently expressive such that they can be used to solve the targeted application (embedding) problem. We then sample from the input space X using a statistically grounded sampling method (e.g., Latin Hypercube sampling). For output points which map into the (pre-defined) high information content region Y * , we numerically approximate the associated gradient information (e.g., via (24)), and we apply the directed walk method (e.g., via (25)) to enrich the training dataset. Gradients associated with all data are then computed in order to regularize the training process via (12).
Before training, the dataset (with the directed walk points) is split into training, validation, and testing subsets and the input and output datasets are standardized. Then the training is conducted using the validation dataset as a monitoring of the generalization error. For the "loop-closing" processes the training is halted at an intermediate stage and one of the two subsequent resampling procedures is initiated. First, a large set of new input samples are chosen and rapidly evaluated with the NN. Points which map intoŶ * are saved and evaluated by the ground truth mapping (2); we refer to this process as Neural Network Informed (NI) sampling. Second, we reformulate the (ReLU-based) NN as the MILP via (4), (5), (7), (10a)-(10e). The resulting MILP is used to find both adversarial points (i.e., points where small perturbation lead to large prediction changes) and points whose predictions maximally deviate from some expectation (see Section IV-D for examples and details). We refer to this process as Verification Informed (VI) sampling. All model training undergoes a hyper-parameter tuning before assessing the performance on the final test dataset. Lastly, the MILP reformulation is used to embed the now-trustworthy model inside a relevant problem (see Section IV-E for examples and details).

IV. CASE STUDY: THE NORTH SEA WIND POWER HUB
In this section, we present a proposed model of the North Sea Wind Power Hub (NSWPH). This model represents a realistic example of an emerging converter-dominated system that must integrate securely into the legacy bulk power grid. Converter-dominated systems have many more degrees of freedom compared to conventional power systems, and they also require several additional control loops to maintain a  Figure 3. The North Sea Wind Power Hub. The offshore AC hub connects to five wind turbines, and a local synchronous condenser (SC) that offers inertia and voltage support. Hub voltage source converters (VSCs) supply DC power to onshore (mainland) power systems via HVDC sub-sea cables. The associated dynamical system model contains 164 differential states. stable operation. Tuning all necessary control parameters to ensure stability and achieve a good performance becomes an increasingly complex task; it requires advanced optimization tools to explicitly consider the stability constraints of the system. This has so far been very difficult to achieve with conventional optimization tools. This paper proposes a machine learning framework to overcome this otherwise intractable optimization problem. The objective is to optimally tune the control parameters and determine the power dispatch of the NSWPH system under small-signal stability and N-1 security constraints. Figure 3 depicts a circuit diagram of a proposed North Sea Wind Power Hub (NSWPH), with 5 wind farms, a synchronous condenser, and two HVDC Voltage Source Converters (VSCs) that connect the offshore AC Hub to the mainland. In modeling this system, we focus exclusively on the dynamics of the offshore AC grid (Hub) that includes the synchronous connection between the wind turbine VSCs, the hub VSCs and the synchronous condenser. The technical details associated with the NSWPH system and VSCs are reviewed in Appendix A and in [1], [78].
When planning the dispatch of this system, operators are assumed to have control over the active and reactive power set points of the turbines (P ref and Q ref , respectively). In future power systems, operators may also depend on the wind hubs to provide both primary frequency and primary voltage control support services. Accordingly, we assume operators have the capacity to adjust the droop parameters which determine the system's participation in both primary frequency and primary voltage control (K p,f and K v , respectively). These parameters are depicted inside the active and reactive power control loops of the VSCs in Fig. 12.
When dispatching the power set points and control regulation, operators must ensure that local system damping ratios ζ i remain above some minimally acceptable threshold (e.g., 3%). In this paper, we assume this margin has to be maintained across all credible N-1 contingencies (i.e., the loss of any one turbine). The mapping between these parameters and the minimum system damping ratio, however, is intractably dense (see Appendix A), since it involves solving a series of N-1 power flow problems, initializing state variables, linearizing the network, and then computing the damping ratios of this linearized model. To help overcome imposing this computational burden in optimal control problems, we develop a NN mapping between the tunable control parameters and the minimum damping ratio of the NSWPH across all N-1 contingencies: whereζ c is the estimated minimum damping across all N-1 contingencies and x,x represent the upper and lower bounds of the set points and parameters. In the following subsections, we present the methods we use to sequentially (i) perform data collection and enrichment, (ii) pre-process the data, (iii) regularize the training procedure, (iv) verify model behavior and improve data sampling, and (v) embed the associated model inside an optimization problem.

A. Dataset creation and enrichment
To collect data for training the NN model of (13), we began by defining a 4-dimensional hypercube and K v across their lower (l) and upper (u) limits. After sampling H on each dimension a total of N times, we mathematically computed the minimum damping ratio ζ c across all N-1 contingencies associated with each operating point x. We tested three types of sampling techniques: • equally spaced grid (i.e., grid search); • uniform sampling in each dimension; • Latin Hypercube (LHC) sampling. For each routine, we sampled 5 4 , 6 4 , and 7 4 data points.
While these sampling approaches can evenly sample across the full input domain, we are most interested in collecting data around the pre-defined stability margin. Accordingly, once we computed the minimum damping ratio associated with the i th sampled operating point x, we performed a directed walk [11]. This directed walk iteratively progressed an initial operating point x (0) to points x (1) , x (2) , etc., each of which mapped to an eigenmode which was incrementally closer to the stability margin ζ c = 3%. Notably, this process was only implemented if the minimum damping ratio of the initial operating point was sufficiently close to the stability margin (e.g., 3% ± 6%). The procedure is depicted in Fig. 4, and the details of its implementation are described in Appendix B. For each point sampled from the hypercube, and for each terminal directed walk point, we saved both the input (K p,f , K v , P ref , Q ref ) and output (ζ c ) data, but we also saved the gradient information ∇ x = ∂x i /∂y i , since it was also used to effectively regularize the training procedure.
We proceeded similarly for the addition of network informed (NI) and verification informed (VI) points according to Algorithm 1 and Algorithm 2, respectively. For the NI points, we interrupted the training procedure after 1 000 epochs and sampled a 1 000 000 possible input data points x in the hypercube. We then evaluated the NN at these points, but retained only predictions within a region close to the stability boundary, defined by theζ c ± ∆, where ∆ determines the adjustable width of the region. We then sampled N samples = 200 additional points from this remaining set. These points were then evaluated using the physical model and ultimately added to the training dataset.
For the VI points, we identified regions of consistent classification starting from each of the 16 corner points of the hypercube H. We then sampled a large number of points in the hypercube, and we omitted those input points which fell in the regions that were apriori verified stable or unstable. From the remaining points, we again sampled N samples = 200 points which we then evaluated. See Section IV-D for the exact implementation of VI sampling.

B. Data pre-processing
As part of the data handling, we classified the points into five classes for the re-sampling and, later, model assessment. These classes are stable, marginally stable, marginal, marginally unstable and unstable, and they are defined relative to the critical damping ratio as shown in Table I. Alongside the definition, we display the distribution of the classes in the test dataset. We used a very large test dataset (N = 21 4 = 194 481) equally spaced in the normalized hypercube. This case may not represent the usual practical setup, as it is rarely feasible to generate such a comprehensive dataset that can be considered as the ground truth. However, we opted for it to show the performance of each model as objectively as possible. Two important notes on this implementation: 1) The test set was only used for "final" evaluation and was not part of any feedback loop or model selection.
2) The assumption of having the same underlying statistical distribution does not exactly hold for the datasets generated using the LHC and uniform sampling methods. Therefore, this exhaustive test dataset favors the dataset also sampled on a grid; this should be kept in mind when interpreting the final results.  For the training/validation split, we opted for an 80/20 split of the datasets, and we fixed how the dataset was split using a constant "seed" for the utilized random number generator, such that all models performed under identical conditions. The sole data manipulation we conducted in the pre-processing stage was to standardize both the input and the output data. For this, we computed the mean and standard deviation of each input and output in the training dataset, and we applied the same standarization values also to the testing and validation datasets.

C. Model training and gradient-based training regularization
The objective, or loss, functions of the training problem are defined as where N stands for the number of data points in the respective dataset. Equation (16) All models were trained with the Adam optimizer [79] with a decaying learning rate l that depends on the epoch n via where l 0 represents the initial learning rate and γ the learning rate decay. Each NN was trained for in total 3'000 epochs, but, in the spirit of "early stopping", the weights and biases that achieved the lowest loss on the validation set were used for model selection purposes. For all model forms, we ran a hyper-parameter exploration using the values shown in Table II. This was performed by conducting a grid search over the entire parameter space; the best hyper-parameter set was selected by using the validation set. Each combination was run ten times with different weight initialization seeds, and the hyper-parameter set that achieved the best mean of the validation losses on the different initializations was selected. The resulting sets of hyper-parameters are summarized in Table III in the appendix. For the final model comparison, each model was run with 100 different initialization using the identified best hyper-parameters.

D. Model assessment and verification formulation
For the statistical model assessment, we utilized the previously introduced large test set (N = 21 4 ). Besides the total loss across the test set, i.e., (15), we also computed the test loss across each of the five classifications C i in Table I by substituting N in (15) only with the N i points that belong to the respective classification.
The verification problem (for VI sampling) sought to find the smallest distance ∈ [0, 1] in the normalized hypercube for which the classification threshold (see Table I) is passed. If a given corner point x 0 was stable, then the threshold was set to 3% + ∆, and if x 0 was unstable, then it was set to 3% − ∆. For marginal and marginally un/stable points, the verification problem was skipped; Fig. 5 depicts the procedure in twodimensional space. The mathematical program associated with a stable corner point is wherex indicates component-wise normalization (x 1 = x 1 /x 1,max ). If the corner point is originally unstable, then (19b) changes toζ c ≥ 3% − ∆. As a comparison, we also ran a data-driven version of (19): instead of using the NN's MILP formulation, we used a sampled collection of points. By approximating in that sampled set of (x,ŷ x ) we could compare the effectiveness of the statistical and continuous evaluation techniques.

E. Embedding problem
By employing the MILP reformulation from (10a)-(10e), the learned NN mapping (13) can be directly embedded inside of a variety of otherwise intractable optimization problems. One such problem which motivates this work is the cooptimization of droop parameters in the context of a regional planning study (e.g., as in [80]). Such planning studies often necessarily utilize overly simplistic linear models of frequency dynamics. Through NN embedding, however, linear models can be replaced by piecewise linear NN-based models which have the capacity to be far more accurate [4].
In the NSWPH case, for any targeted nominal power injection For each of the four corner points (x (c1) through x (c4) ), the MILP (19) attempts to minimize while the infinity norm (green box) touchesζc ± ∆, which is the marginal stability bound predicted by the NN. The MILP thus generates "in/stability certificates" in these regions, guaranteeing how the NN will classify all associated internal points. The certificates are then used to filter out future samples. No certificates are generated at x (c2) , since the stability boundary overlaps closely with the corner point itself. For reference, we also show * 1 in the lower box, which is the ground truth value of 1 .
regions of droop parameters around a given parameter setting K p,f,0 , K v,0 which are predicted to yield definitely-stable system behavior: The solutions of these embedding problems are illustrated in the subplots of Fig. 10.

F. Notes on the implementation
The entire framework, which is provided on our Github [9], is implemented in Python 3.8 using standard libraries. All parts related to the training of the model utilize PyTorch [81]. The optimization problems are formulated in Pyomo [82], [83] and solved using Gurobi [84]. The logging and tracking of all training runs, in particular for the hyperparameter tuning, is conducted with WandB [85]. The training and optimization were performed on the machines in the high performance computing (HPC) cluster at the Technical University of Denmark (DTU) [86] that are configured as nodes with 2xIntel Xeon Processor 2650v4 (12 core, 2.20GHz) and 256 GB memory.

V. RESULTS
The following paragraphs show the previously elaborated steps in greater detail by illustrating important aspects of the procedure. It is important to realize that we encountered competing objectives when training the NN: on one hand, we aimed to achieve an accurate predictor across the full input domain in order to use the NN as a representative function in an eventual embedding problem. On the other hand, we tried to capture the region around the critical damping ratio as closely as possible. With these competing objective in mind, we used data associated with the full input space as the decision criterion for all model selection decisions, but by design, the methods that sampled additional training points targeted data close to the stability boundary.

A. Problem characteristics -class imbalance
Before diving into the direct comparison of the methods, we take a brief look at the characteristics of the problem, namely the mapping f from the input to the output domain. The fundamental problem in the dataset creation becomes apparent in Fig. 6, which illustrates the relationship between the input parameter K p,f , the damping ratio ζ, and the number of data points collected in each region. With simple sampling methods (which are used to generate the upper row in Fig. 6), the vast majority of data points lie well within the stable region; this also means that very few data points fall in or close to the marginal region. In contrast, a very high percentage of the DW dataset (lower left) hits the marginal stability region, as it was designed to do. Similarly, the NI resampling (∆ = 3%) creates a much more targeted sample distribution. Lastly, the VI resampling (∆ = 0.25%) is not as successful in creating a more balanced distribution, mainly due to the fact that it only can omit about 40% of the input domain. Table IV in the Appendix reports a more detailed overview of this class imbalance, and it shows how the resampling distribution varies across different training runs associated with the NI and VI approaches.
Understanding these problem specific dataset characteristics is a crucial step for the subsequent performance assessment and analysis. This also highlights why the data pre-processing step should also include thorough data exploration; simple tools, such as the regional classification of Fig. 6 or correlation analysis, can already provide vital insights.

B. Comparison of the model performance
With these considerations, we move to the central results of this work: the model assessment in order to compare the effects of the previously presented approaches. the different variations of the models for the two primary metrics of interest: the overall mean squared error (loss), and the loss in the marginal stability regions that we target to improve. The Box-Whisker plots stem from NN training repetitions across different initializations, subsequently called "runs", of the NN's weights and biases. Across the different sampled base datasets, it becomes clear that simply adding more data points improves the test accuracy (the loss on the test set) in terms of the mean, but also for the standard deviation, across the different runs. This is simply the effect of better model generalization and less dependence on a good training run by virtue of a "lucky" initialization. The addition of the physics regularization (PR) shows a similar effect for the datasets stemming from LHC and uniform sampling. This is a clear sign for an effective regularizer -in contrast, the PR for the Grid dataset is much less powerful. This specific effectiveness, which depends on the sample distribution, may stem from the following consideration: the gradient regularization helps the NN interpolate better in between data points. However, for the Grid dataset, the data points are already implicitly carrying this gradient information, albeit not very accurately. Hence, the interpolation to points that lie not on the lattice is not very much aided by the additional regularizer, i.e., the additional information content is not as high as for the other sampling strategies.
We have to note that these results are not sufficient to draw conclusions about the best sampling approach, as only a single dataset was created for the LHC and Uniform approach. For conclusions on the sampling efficiency, these tests would need to be repeated several times with different realizations of the sampled datasets. Nonetheless, it is striking that both the NI and VI approach yield more consistent results, i.e., smaller standard deviation, for LHC and Uniform sampling than for the grid dataset. The cause probably lies in the much less accurate predictions of some training runs at the time when the additional points are sampled. Hence, the additional points will not be as close to the target region, which leads to lower performance improvement. In this regard, we can also observe that the NI approach has an advantage over the VI approach across all sampling strategies when measured by the test loss in the marginal region. This is not surprising, as we showed in Fig. 6: the NI approach is much more effective at sampling in the marginal region. The VI approach can reduce the size of the input domain, but mostly only by a 40-50% measure in terms of the volume of the unit hypercube. Table IV in the Appendix supports this observation, by analyzing the share of the (additional) samples of each class. Unsurprisingly, the DWs yield the best performance among the "resampling" strategies when considering the loss in the marginal region. This is particularly noticeable for the Grid dataset, as the DW outperforms both the NI and VI approaches.
Ultimately, these analyses show that by introducing additional elements into the ML framework, in particular feedback loops, we can create well performing models without the need for greatly increased sampling efforts. To judge the efficiency of these methods, we consider how well the simple sampling strategies perform given an increased number of points (approximately double and four times for datasets with 6 4 and 7 4 sampled points). Whereas the performance improvements for DW, NI, and VI on the overall test loss can be attributed, to some extent, to the additional number of data points (usually around 200), the test loss in the marginal region clearly shows that a target sampling of these ca. 200 points can lead to performances that would need significantly more additional data points with conventional sampling methods. Furthermore, if an effective regularizer can be provided, such performance boosts can be obtained without additional data. Table V in the Appendix presents all the results in a numerical form.

C. Hyper-parameter tuning and model selection
The previously presented results represent only the model assessment step; the model learning and, in particular, the model selection steps are often given much less attention. Therefore, we want to stress the importance of a thorough hyper-parameter tuning and model selection procedure. The presented framework and its implementation allow for easily running a large number of experiments with the sole restriction being the computational burden. The HPC cluster we used at DTU allowed us to systematically test all combinations of hyper-parameters listed in Table II. After running these combinations with 10 different initialization seeds each, we used the mean of the loss across the respective validation set as the model selection criterion, which led to the selection as presented in Table III  relevant. For the model to generalize well, model selection must not happen based on the test set that is used for the performance assessment -otherwise, the assessment will be biased. Having distributionally invariant subsets of the datasets will allow us to choose well performing models nonetheless. In Fig. 8, we show what the performance assessment distributions that the different hyper-parameter combinations would have yielded (boxplots), together with the hyper-parameter setting that was chosen based on the above mentioned model selection criterion (red dot), i.e., the loss on the validation set. We can observe that in all cases the overall test loss is close to the best possible model, and similarly for the test loss in the marginal region. However, there are a few instances where there could have been a better hyper-parameter selection, e.g., for the Grid VI model. This insight, however, is only an analysis in hindsight, and it must not influence the model selection process. On the other hand, we can also see how a poor hyper-parameter selection could have led to much worse model performance and hence completely cancel out the performance benefits of the various methods proposed in this paper and invalidate all previously presented results.
It should be noted that the extent of this hyper-parameter study is usually beyond the scope of a practical implementation. However, for a reliable performance comparison, in particular if conclusions about the superiority of a method are to be drawn, a sufficient exploration of the utilized hyperparameters and a rigorous model selection procedure must be performed. Working in a fixed framework, as presented, can greatly reduce the effort to do so and at the same time place safeguards against (even unconsciously) biasing the results.  Table II). The red dot shows the selected model based on the validation loss of the particular dataset (see Table III). The fact that the selected hyper-parameters do not always achieve the best performance when evaluated on the test set underlines the importance of the dataset splitting. Same coloring as in Fig. 7.

D. Statistical vs. continuous model assessment
A central element in the framework is the reformulation of the NN as a MILP. We use it for the VI resampling, and for the embedding problem, and it can be used also in performance assessment. To show the benefit that we can draw from the ACCEPTED FOR PRESENTATION IN 11TH BULK POWER SYSTEMS DYNAMICS AND CONTROL SYMPOSIUM (IREP 2022), JULY 25-30, 2022, BANFF, CANADA 14 reformulation, we look at how a statistical evaluation (i.e., simply passing data points through the NN to assess stability regions) compares to the MILP verification (i.e., using the MILP reformulation to rigorously assess stability regions).
As an example, we look at the calculation of a certified region starting from two corner points, namely x 1 = [0.0, −0.5, 0.0, 0.0] and x 2 = [0.0, −0.5, 0.0, 50.0]; the former should result in a stable region and the latter in an unstable one. We run the verification according to (19) and subsequently compare the result to the size of the region that a statistical evaluation would yield for increasing the number of tested points. Figure 9 shows the results, here plotted as the distribution across 100 different samplings. It is intuitive that the gap between the verified behavior ( verified ) and its empirical estimation ( statistical ) closes with the increase in the number of sampled points. However, the plot should also make clear that the empirical value can change drastically with the addition of only a single point. It so happened for the blue curve in the lower panel of Fig. 9: around the 1100th point, the gap abruptly diminished to then stay constant until the millionth point. In contrast, the gap for the orange sample only reaches this level after many hundreds of thousands of points are used. Furthermore, a good estimation for one gap does not translate into a good estimation for any other gap, as we can observe when comparing the "good" estimation of the orange sample in Fig. 9a with the "bad" one in Fig. 9b.
This example illustrates that relying on empirical evaluations of the behavior of neural networks requires a large number of points to become statistically trustworthy. And even then, there is still a chance for very small regions that yield a contradicting prediction -only by verifying a NN with optimization tools can we be certain about its definite behavior. In case such a contradicting prediction occurs, we can then furthermore compare it to the ground truth and identify it either as an adversarial example or confirm the seemingly contradicting prediction.

E. Embedding possibilities
As a last topic of the results we come back to the motivation of this work: training a NN as a surrogate function which we then can embed into other problem such as described in Section IV-E. Figure 10 presents two illustrative results for each of the posed problems. Figures 10a and 10b show the verified regions of control parameters that lead to stable behavior for a fixed power injection P ref , Q ref . In the inverted setting, where droop parameters K p,f , K v are fixed, Figs. 10c and 10d identify admissible power injections which could be included, for example, in an OPF setting.
Following up on the previous section of statistical versus verified behavior, the shown regions outlined by the dashed white lines are verified, whereas the contours in the plot stem from NN predictions evaluated on an extremely fine grid (akin to a form of statistical evaluation). Either of the options could be used; however, as the input problem's dimensionality increases, the latter becomes harder and harder to realize.

VI. DISCUSSION
We view the proposed framework, along with the considerations presented throughout this paper, as important steps towards bringing Power Systems and Machine Learning (ML) closer together and enabling the use of ML in safety-critical applications. Borrowing from the ML community, one could find plenty of generic workflows and implemented frameworks and [87] reviews a number of case studies -however, for a successful adoption of ML methods to power systems, we need to tailor these approaches to our specific needs. At the same time, the sheer number of specific problems in the power systems domain makes generic solutions, as well as communitywide benchmarks, difficult to implement, especially since the primary goal of using ML in power systems is still driven by specific applications. We should, nonetheless, strive towards clarity, comparability, and transparency in our efforts in ML, because the alternative (e.g., unclear assumptions, improper dataset handling, or badly tuned algorithms) will certainly not strengthen efforts to adopt ML in safety-critical applications. To establish such practices, we believe that three pillars will be crucial: i) code, ii) data, and iii) communication. In all of these aspects, a framework, such as the one proposed in this paper, can be helpful, and the following paragraphs will briefly discuss how. research, the produced code should be made publicly available; this not only helps demonstrate the strengths of the proposed algorithm (which is very difficult to prove otherwise), but it can also serves as a benchmark for other researchers to compare with. Second, it should not be underestimated how much a framework can help structure the code and thereby make it more adaptable. Ideally, each of modules shown in Fig. 2 is treated as a stand-alone piece of code. This not only allows one to interchange methods more easily, e.g., the sampling method for the dataset creation, but it also allows one to be more aware of the information flow between modules. This means that each of the arrows in Fig. 2 is associated with passing on some information; this can range from datasets to parameters to performance assessments. By being more conscious about which information is being passed on and by trying to reduce these exchanges to the necessary minimum, we become less prone to introducing biases in the process and it gets easier to check for such undesired information flows. The classic example is that the test dataset should not enter the model learning but only the model assessment. And especially if we introduce more loops in the framework, it must be ensured that, e.g., the VI or NI approach do not make use of the test dataset. A side-effect of modularizing code is that the workflow becomes more and more generic which immensely simplifies conducting hyper-parameter studies and the comparison of methods. In the provided code base [9], a small configuration file is used to setup and run all the different setups that are reported. Thereby, the risk of errors greatly reduces and the repeatability of experiments in case of changes in a module is significantly simplified.
ii) Datasets: It should be apparent by now that the used datasets are of utmost importance. Providing access to the models that have been used for the dataset creation (or at least, the resulting datasets) will become a necessity. In case proprietary reasons prevent a release, mock datasets that resemble the used datasets should be tested and provided. Furthermore, it is advised to use so-called seeds in all processes associated with random numbers, e.g., for the train/validation/test split. It ensures the reproducability of the results, and it reminds us to consider how general the drawn conclusions can be. In the presented results, the LHC and uniform sampling approaches have only been performed once, as tests with multiple sample realizations would have increased the computational burden even further. This design choice, which can be challenged, needs to be accounted for when comparing the results across the different sampling methods. A conclusion that would claim the superiority of, e.g., uniform over LHC sampling, is simply not trustworthy based on our results, as this conclusion might have been an artifact of the randomly chosen datasets.
iii) Communication: One of the biggest challenges is to efficiently communicate at the intersection of ML and power systems. It begins with terminology: some of the terms are community specific, some overlap, and some are over-loaded, meaning that the same word can take multiple meanings. ACCEPTED FOR PRESENTATION IN 11TH BULK POWER SYSTEMS DYNAMICS AND CONTROL SYMPOSIUM (IREP 2022), JULY 25-30, 2022, BANFF, CANADA 16 A simple example is the word 'model': we tried to avoid confusion by using NSWPH model or physical model for the power system model and ML model for the learned neural network, but still it is not always easy to be clear. Furthermore, with all additional terminology, there is a degree of ambiguity and it will take time to agree upon how the power systems community will use new terms introduced from the machine learning jargon. For such a setting, a framework as proposed can be helpful as a starting point. Beyond simple terminology, building trustworthy ML also requires a common understanding of what needs to be presented in a ML related research article. By requiring, e.g., an overview of the hyper-parameter tuning or basic statistical properties of the dataset, we might be able to establish good practices and raise awareness of the associated issues among the entire community.
In terms of technical aspects of the presented work, we want to briefly discuss two aspects that have not been addressed in the Results Sections: The overall computational effort and implications for practical work.
On the aspect of the computational effort, it is important to distinguish between the entire workflow and individual steps involved. From a high-level perspective, there is a trade-off to be made between model performance and computational effort. The simplest example of this is for instance that with the possibility to generate additional data points, it should be possible to create a better performing model, but at the cost of computing the additional data points. Therefore, an honest comparison of methods pays tribute to this trade-off. This can, for example, be achieved by reporting the sensitivity of a method to the dataset size, i.e., by training with varied dataset sizes, as we did. If a method, such as the NI approach, can reach a similar performance improvement as a simple doubling of the dataset size would while only requiring a third of the necessary additional simulations, we can consider the method a true improvement. This comparison was rather straightforward as it did not involve measuring how long it takes to simulate a data point. As soon as we enter comparisons between methods that involve computationally intense tasks, the actual implementation becomes a source of uncertainty that is hard to control. For example, the dataset creation can be sped up significantly by performing all calculations in Ccompiled code. On the other hand, the VI approach is strongly dependent on the NN size but also its sparsity. The resulting assessment of whether a larger dataset or the VI approach is computationally more efficient becomes arbitrary as long as these comparisons are not performed with open-source stateof-the-art implementations. To alleviate this complication, one can instead revert to qualitatively analyzing the share of each module or step in the process to identify the biggest burdens and how the proposed methods affect these. In the presented study, the dataset creation clearly dominates the burden of the training (also including the physics-regularisation). For the verification, the picture is mixed, as it quickly requires the biggest computational effort for the unfavorable cases of large NNs. But as mentioned previously, with an improved setup, i.e., a sparser NN, this relation could easily be proven wrong. Accordingly, we did not strive for minimal computation time, and, hence, did not place emphasis on analyzing it in detail. Another, more formal alternative can be to focus on analyzing algorithms in terms of their algorithmic complexity such that their scalability can be judged.
For a practical implementation of the framework, the extent of the hyper-parameter study can be significantly reduced. Instead of a full grid search, which we conducted for a complete picture, it is much more practical to first run a few tests with random hyper-parameters drawn from a larger range before moving to a reduced range with a more thorough grid search. In the study, the low learning rates and large learning rate decays could have for example been discarded early on. A slightly advanced method is to employ Bayesian hyper-parameter tuning. These comments, though, shall not encourage one to take short-cuts on hyper-parameter tuning, but instead conduct it in an efficient manner. A final note on how the suggested methods might be used: From a methodological point of view, the physics-regularization and the NI resampling approach appear to be the most promising at the moment. The former appeals due to its simplicity and avoidance of the (expensive) creation of additional data points, the restriction being that the required values, here derivatives, are available. On the other hand, the NI resampling approach stands emblematic for the idea of "closing the loop" as it combines the cheap evaluation and selection of candidates via the NN and a targeted dataset enrichment strategy. While not as targeted as the directed walks, it achieves a similar effect without the majority of its cost. Regarding the VI approach, improving the size of the rejected areas, e.g., by including points other than the corner points, will be crucial to use it for general resampling purposes. The true potential of verificationbased approaches lies in the structured search for adversarial examples and the possibility to provide regions of verified behavior for the embedding problems; these are elements that a purely statistical assessment cannot provide.

VII. CONCLUSION
In this paper, we identified key challenges associated with building high-fidelity ML models for power system applications. We then integrated state-of-the-art solutions to these various challenges into a single coherent framework. The proposed workflow at the heart of this framework utilizes tools such as physics regularization, hyper-parameter sweeping and loop-closing in order to efficiently generate ML models which are both trustworthy and generalizable. We tested the proposed framework on a realistic model of the North Sea Wind Power Hub. Through extensive analysis, we quantitatively demonstrated the benefits that can be provided by the various elements of the proposed framework. We also illustrated how derived ML models can be used to generate useful embedding (e.g., (20)) and interpretable mappings (e.g., Within the next decade, governments and grid operators around the North Sea plan to install offshore wind power hubs with never-before-seen capacities. The study on the "North Sea Wind Power Hub" (NSWPH) in [1] details one possible configuration of this system. The AC dynamics of this system will be dominated by the complex dynamics of Voltage Source Converters (VSCs). Figure 11 shows the inner control loops on the AC side of a typical VSC [1], [78]. The active and reactive power controllers associated with this VSC are depicted in greater detail in Fig. 12. In particular, these power controllers explicitly depend on droop parameters which determine the system's participation in both primary frequency and primary voltage control (K p,f and K v , respectively). The influence of these parameters inside of the active and reactive power VSC control loops is depicted in Fig. 12.
As shown, the control loops also depend on active and reactive power reference set points as input parameters. In the optimal dispatch problem, we assume grid operators have simultaneous control over P ref , Q ref , K p,f and K v . Notably, K p,f and K v only affect the dynamics of the NSWPH system (rather than its steady state behavior). The active and reactive power set points, however, affect the power flow initialization of the entire circuit, thus implicitly altering the system's linearized dynamics.
When dispatching the power set points and the quantity of control regulation provided by this system, operators must ensure that local system damping ratios remain above some minimally acceptable threshold (e.g., 3%). The damping ratio ζ i associated with the i th complex eigenvalue λ i = σ i + jω i is given by In this paper, we assume this margin has to be maintained 2 across all credible N-1 contingencies (i.e., the loss of any one turbine). The mapping between the input parameters P ref , Q ref , K p,f and K v and the minimum system damping ratio, however, is intractably dense; for any change in their values,

B. Directed Walk Procedure
To perform the directed walk, we follow a procedure similar to the one outlined in [11]. For each operating point, we compute the sensitivity of the n th eigenvalue to the i th element x i of the operating point [11]: where ψ T n and φ are the left and right eigenvectors, normalized such that ψ T n φ n = 1, ∀n, and A x is the state matrix of the NSWPH in Fig. 3, linearized around operating point x. Using the real and imaginary sensitivities in (23), damping ratio sensitivity is given as ∂ζ n ∂x i = ω n σ n ∂ωn ∂xi − ω n ∂σn ∂xi The concatenated gradient information ∇ x allows us to push the operating point x towards the stability boundary via steepest descent: where the magnitude of the stepsize α(x (i) ) is adaptively chosen based on the distance to the stability margin, and its sign (+/−) is chosen based on which side of the stability margin the operating point x (i) is on.
In this work, we utilize a descent procedure which holds control parameters (K p,f , K v ) constant while the power set points (P ref , Q ref ) are iteratively updated. Once the stability margin has been reached, the initial operating point is returned to, and the procedure is reversed (power is held constant while control parameters walk).   Table IV   SIZES OF THE STABILITY CLASSES IN THE USED DATASETS TO ILLUSTRATE THE CLASS IMBALANCE. FOR THE VERSIONS WITH RESAMPLING (+ DW, + NI, AND + VI) ONLY  THE ADDITIONALLY SAMPLED POINTS ARE REPORTED. DUE TO THE DEPENDENCY ON THE INITIAL TRAINING, THE ADDED POINTS FOR THE DATASET + NI, + VI VARY FOR   DIFFERENT TRAINING RUNS, HENCE THE MEAN AND STANDARD DEVIATION OF WHICH SHARE BELONGS TO EACH CLASS ARE