Prosumer Response Estimation Using SINDyc in Conjunction with Markov-Chain Monte-Carlo Sampling

: Smart grid operation schemes can integrate prosumers by offering economic rewards in exchange for the desired response. In order to activate prosumers appropriately, such operation schemes require models of the dynamic uncertain price-response relationships. In this study, we combine the system identiﬁcation of nonlinear dynamics with control (SINDyc) algorithm with Bayesian inference techniques based on Markov-chain Monte-Carlo sampling. We demonstrate this combination of two algorithms on an exemplary system in order to obtain parsimonious models alongside parameter uncertainty estimates. The precision of the identiﬁed models depends on the identiﬁcation experiment and the parameterization of the algorithms. Such models may characterize the prosumer response and its uncertainty, thereby facilitating the integration of such entities into smart grid operation schemes.


Introduction
Prosumer response (PR) activation, similarly as DR, is considered a key ingredient in a smart energy system [1,2]. A prosumer is a unit within the power system that can act both as consumer and producer. Examples are EV with V2G functionality. Such cars can extract and inject power from and to the power system [3]. Inclusion of prosumers into the system operation consequently leverages flexibility potentials, which may facilitate the integration of higher shares of RES, thereby contributing to a more sustainable grid operation [4,5]. Furthermore, this may lead to a reduction of the cost of system operation [6]. In order to coordinate prosumers during real-time operation in an optimized manner, a control scheme requires knowledge of the PR dynamics. Such a control scheme can act pro-actively, contrary to naive flexibility schemes; see for example [7].
System identification (Sys-ID) techniques are one central building block for achieving long-term reliable control; see e.g., [8][9][10]. By treating the system as a black-box, Sys-ID enables the control of systems without explicitly modeling the underlying physical properties. The SINDyc algorithm [11] is a recent addition to the Sys-ID field. It builds on sparsity-promoting optimization techniques such as the LASSO algorithm [12]. A response model derived using SINDyc is sparse in the model coefficients. SINDyc therefore aims at an accurate model with a low number of active terms selected from a candidate model structure. A model with such properties is referred to in the literature as parsimonious [13]. Furthermore, SINDyc models have been shown to perform well when facing scarce availability of data [14].
This data-driven approach bears potential for smart grid applications where, for example, privacy restrictions apply [15]. The application of Indirect Control (ICo) is one approach to the integration

A Sparse Nonlinear System Identification Algorithm
The work in [11] formulated the so-called SINDyc algorithm for the identification of sparse nonlinear models subject to forcing, based on the SINDy algorithm [13] for the identification of autonomous nonlinear systems.
In general, we can state the nonlinear price-sensitive dynamical system as: f is consequently governed by free and forced dynamics. We may reformulate it as: Hereby, Ξ are sparse model coefficients. Θ T is the model structure, that is, model terms corresponding to Ξ. X is the Sys-ID data, that is system observations recorded from an Sys-ID experiment. As we consider a system excited by the price-signal p, we split X into the system state x and the price signal p: x is hereby approximated using the variation over Sys-ID data x. In the simplest form, a one-step shifted version of the input-output observations x subtracted from the original version yields the approximations of the dynamics. The related dynamic mode decomposition (DMD) algorithm [33,34] uses a similar approach, however for the identification of linear dynamical system models. Θ T (x, p) denotes the model structural terms including the state x, the input p, and potentially cross-terms of both x and p.
The choice of the model structure is one important design choice [8,9]. A straightforward assumption is to assume Θ T (x, p) to resemble the power series up to a chosen degree. The work in [9] outlined the drawbacks of this model structure resulting from the properties associated with polynomials: • Structure selection is computationally demanding, especially for high dimensional problems.

•
The extrapolation capabilities of the power series are sub-optimal. • Polynomial models suffer heavily from the curse of dimensionality.
Positive properties include [9]: • The capability to approximate a broad group of target problems; • Low sensitivity to noise; • Global explanatory capabilities.
Referring to discussions on this type of model in [9], we remark that we chose this model type for the purpose of demonstrating the application of the SINDyc algorithm in conjunction with MCMC. Other applications may require a different model type. As recommended in [9] and in order to limit the aforementioned drawbacks, we only considered polynomials up to third order. We here use the model structure given in [11], a power series including cross-terms.
Ξ is obtained using SINDyc based on the sequential thresholded least-squares algorithm proposed in [13]. See Algorithm 1. SINDyc_params are parameters passed to the SINDyc algorithm. As outlined in [11], we have to choose the regularization weight α in order to obtain a sparse model whilst retaining model accuracy. We here perform a naive sweep over a set of candidate weightsᾱ as suggested in [11] while evaluating the sparsity of Ξ. We evaluate the sparsity based on the count of nonzeros and the count of values in Ξ and compare it to a chosen sparsity threshold sparsity_threshold. Furthermore, we evaluate the model using a model evaluation function denoted evaluation_function. The model evaluation function should as a fundamental task examine the stability of the model for a set of operating points and excitation signals. It may include a set of additional evaluation tasks.

Probabilistic Model
We may pose (4) as a probabilistic model and state it as: m denotes the prior, which includes the model coefficients Ξ obtained using SINDyc in Algorithm 1: Following the Bayesian principle, we can state the prior in a flexible manner based on available information. Parts of the prior may be undefined. Such a lack of information becomes part of the overall uncertainty in the model. We include weakly informative priors for these parts as recommended in [35].

Probabilistic Model Inference
The inference process of the probabilistic model (5) is formulated in Algorithm 2. XI is a list in which we aggregate models inferred using Algorithm 1.X is a list of individual Sys-ID experiments. For each datum X inX, we call Algorithm 1 and obtain a corresponding candidate model Ξ. select_Xistar selects the MCMC candidate model Ξ based on the collection of candidate models XI. MCMC_function then performs MCMC using the candidate model Ξ . MCMC_params are parameters for the MCMC_function. Algorithm 2 returnsΞ, the posterior PDF of the model coefficients.
Algorithm 2: Probabilistic model inference: using a candidate model Ξ derived using SINDyc in Algorithm 1, MCMC is performed on the system observationsX.

Excitation Model
We used the software package Stan [31] for performing MCMC. Stan requires an ODE with modeled forcing for the inference of the dynamics subject to forcing. We therefore augmented the system in (4) with a forcing model that approximated the excitation signal. We restricted the excitation model to a third order polynomial as recommended in [9].

Results
Consider a system of two prosumers. The first prosumer is governed by nonlinear dynamics; the second prosumer is governed by linear dynamics: We observe the system response through the measurements y subject to white noise v as: ω is the noise in the dynamics as introduced in (5). The scalar p is the price-signal sent to the prosumers. We parameterized the system as given in Table A1. We chose a sampling rate T s = 1s. We considered two clusters of data c = {c 0 , c 1 }, where c 0 contains n obs,c 0 = 50 independent system observations series and c 1 contains n obs,c 1 = 5 independent system observation series. While this is a simple model, it should suffice to outline the modeling approaches described in the following.
As for Sys-ID in general, the choice of the excitation signal is fundamental for the quality of the system approximation [8,9]. The excitation signal should hereby correspond to the magnitude and frequency range in which we aim to use the model [9]. Whether the excitation signal is adequate to extract sufficient information is to be checked in relation to the considered system and its operating condition. The sparsity structure is one criterion to determine the success of the Sys-ID experiment.
Here, we chose a double-sinusoidal excitation signal p e applied on top of an assumed constant signal p 0 = 0.5. The constant signal excites the balanced system throughout a burn-in period, such that the system settles on a new approximate equilibrium prior to the start of the system identification period: The criterion chosen to identify the convergence of the system to this equilibrium must be of an approximate nature and in correspondence to the system uncertainty.

Polynomial Prediction Model
We aimed at a low order model as the simplest candidate model without drift term, such that the system remains in equilibrium when unexcited. We chose the candidate model structure as: For the SINDyc algorithm, we choseᾱ as the n α candidate regularization coefficients α equidistant in the sweep bounds b. See the parameters in Table A1 and Algorithm 1.
Identifying models using Algorithm 1, we obtained the coefficient distributions for each cluster c 0 and c 1 illustrated in Figure 1. The uncertainty in the dynamics α m and β led to uncertainty in the magnitudes of the model coefficients.
Model m 0 was correctly associated with nonlinear dynamics, and Model m 0 was correctly associated with linear dynamics. This only held when choosing a well-posed excitation sequence. The convergence of the SINDyc algorithm was assured only when the identified model was sparse in its coefficients. Consequently, the examination of the sparsity of the identified model provided information on the success of the identification. The uncertainty in the parameters was higher for cluster c 1 , in which we had access to only five system observation sequences. The SINDyc model of the first PR approximates the true system for both the identification data and when considering an out-of-sample experiment. See Figure 2. This works in a similar manner for Model m 0 and also for the second cluster c 1 . Notice that the design of the out-of-sample excitation signal, in relation to the Sys-ID excitation signal, is an important step to assess the flexibility of the model in describing the true system. We can visually examine the quality of the fit by comparing the one-step prediction surfaces of both the true system and the deduced model. See

Model Coefficient Distribution Inference Using MCMC
We now aimed to obtain a probabilistic dynamic system model of the first prosumer based on the identified candidate Model m 0 depicted in Figure 1.
At least a subset of the n c 0 observations in cluster c 0 resulted in feasible candidates models, which may inform the prior PDF of the probabilistic model. We may neglect all zero-terms in the candidate model structure in (8) such that the inference through MCMC used only the candidate model as stated in (8). This reduced the computational burden in MCMC. As stated previously, MCMC can handle more complicated prior PDF formulations. The modeler can state more complicated priors based on available system knowledge. For the sake of simplicity, we used Gaussian priors in the exemplary prior PDF stated below.

An Exemplary Prior PDF
We may describe the observation z through the model output y with normally distributed measurement error with variance σ 2 y , corresponding to the measurement equation of the true system stated in (6b): A log-normal distribution for σ 2 y may be reasonable, such that the mean is log(σ 2 y0 ). Assuming a constant measurement noise magnitude, we could hereby choose the logarithm of the variance of the observations in the burn-in period. Notice that for the parameter values chosen here, the Gaussian prior model for σ 2 is appropriate. However, for some prior of σ 2 , it can be appropriate to use the natural conjugated prior, the inverse Gamma distribution. See [36] for examples.
We assumed the prior for the model coefficients Ξ to be normally distributed around Ξ c0 , the coefficients inferred using the SINDyc algorithm:Ξ The lack of information in the formulation of this prior we may express through statements of weakly informative priors; see [35] for examples. Furthermore, we considered the following PPC for the evaluation of the accuracy of the inference:

Probabilistic Model
Stan solved the ODE considered here using a Runge-Kutta-45 method. For MCMC, we chose n iter iterations per chain, n chain chains, n burnin burn-in or warm-up iterations, no thinning, and a seed specified as ς. See the parameter values in Table A1, and consider also [31] in this regard.
As the kernel density estimation bandwidth, we used Scott's rule as given in [37]. The approximated posterior distributionsΞ c0 are depicted in Figure 4a alongside the true parameter distributions and the parameter mean inferred by SINDyc. We can observe that the inferred parametric uncertainty was comparably reduced in cluster c 0 (see Figure 4a) in comparison to cluster c 1 (see Figure 4b). The posterior PDF deviated from the true system's parametric distributions. Hereby, the parametric mean of the SINDyc models provided a higher precision with respect to the true system's parametric mean than the inferred parameters.  Improving the prior optimizes the sampling space for a new MCMC sampling. This can lead to shorter computation time for future sampling of the model. Examining the PPC illustrated in Figure 5a reveals that the model approximated the observed output sub-optimally, yet captured the general trend in the data. By means of random draws from the posterior samples, we could obtain a probabilistic model. Out-of-sample co-simulation of this model alongside the five samples of the true system is depicted in Figure 6.  Quantile range around median

Discussion
Activation of system flexibility through DR and PR schemes such as ICo is increasingly relevant in relation to power system operation. See [39] for an example.
Bayesian approaches are computationally complex [26]. The combination of sparse system identification algorithms such as the SINDyc algorithm and MCMC facilitate the inference of model parameters alongside parameter probability estimates. As shown in the related publications [26][27][28], Bayesian inference techniques may benefit from sparsity, promoting modeling approaches. The focus on a viable candidate model and associated parameter spaces reduces the problem size. Aside from this, parsimonious candidate model can be a desirable goal in the modeling process [13].
Similarly, as described in [27], SINDyc models may provide information on whether nonlinearities are present in the data. This works when the algorithm is provided with informative data and a reasonable library of candidate models. Contrary to [27], we coupled the SINDyc algorithm with MCMC in order to obtain a prior PDF that could reduce the computational burden. Automatized Sys-ID pipelines could benefit from such information.
We designed the excitation sequence in the Sys-ID experiment naively by evaluation of the sparsity of the derived SINDyc models. Due to the fact that such experiments must, at least partly, take place during regular system operation, the modeler may consider active learning techniques. Then, the excitation sequence can be designed such that also operational requirements are accounted for. See [40] for an example. Similarly important as the design of the Sys-ID experiment is the design of out-of-sample simulations for evaluation purposes in order to evaluate the model against the true system.
We chose the SINDyc models based on the evaluation of the combined parametric likelihood. The single candidate model then provided the basis for the formulation of the prior PDF. As depicted in Figure 4, the chosen SINDyc parametric mean values introduced a loss in precision with respect to the true mean values. This was a realistic assumption, due to the fact that the true parametric distributions are unknown in the real setting. Similarly, the inferred parametric posterior PDF were sub-optimal with respect to the true distributions.
Hereby, the parametric mean of the SINDyc models provided a higher precision with respect to the true system's parametric mean than the inferred parameters. In conjunction with the PPC depicted in Figure 5, this may point towards a sub-optimal prior formulation. Furthermore, the polynomial excitation model of third order led to a loss of precision. See Section 2.4 (excitation model).

Materials and Methods
The code used to generate the results presented in this paper can be obtained through the following public repository: https://lab.compute.dtu.dk/freba/sindyc-and-mcmc-framework.

Conclusions
In this paper, we presented a combination of the SINDyc algorithm and MCMC in the context of PR estimation. While SINDyc identifies sparse and potentially nonlinear dynamical system models, MCMC enables the estimation of potentially complex posterior distributions. MCMC can hereby use a sparse system model identified using SINDyc. Due to the fact that sampling is in this manner performed on a constrained and well-posed sampling space, convergence in MCMC benefits from parsimonious models, and the computational load may be reduced. SINDyc may yield such models, provided the appropriate formulation of the system identification experiment and parameterization of the algorithm. Probabilistic dynamical system models enable the application of SMPC, a core ingredient when aiming to activate prosumer dynamics based on informative grounds.
For future improvements, we may replace the polynomial excitation model described in Section 2.4 with an alternative candidate model structure. Alternatively, we may design the Sys-ID experiment such that we can represent the excitation sequence with the low-order polynomial excitation model described herein. In any case, the objective along these lines is to achieve a high accuracy representation of the excitation signal while maintaining a high performance of the sampling process within the MCMC framework. Furthermore, the combination of SINDyc and MCMC should be evaluated on a broad range of modeling problems. When aiming for automatized Sys-ID, robustness and associated issues are to be investigated and potential solutions to be examined.

Conflicts of Interest:
The authors declare no conflict of interest.  Table A1. Scenario parameterization.