Prediction regions based on dissimilarity functions

ISA Transactions 139 (2023) 49–59 Contents lists available at ScienceDirect ISA Transactions journal homepage: www.elsevier.com/locate/isatrans Research article Prediction regions based on dissimilarity functions A.D. Carnerero a,∗, D.R. Ramirez a, S. Lucia b, T. Alamo a aDepartamento de Ingenieria de Sistemas y Automatica, Universidad de Sevilla, Spain bLaboratory of Process Automation Systems, TU Dortmund University, Germany article info Article history: Received 15 July 2022 Received in revised form 13 March 2023 Accepted 31 March 2023 Available online 5 April 2023 Keywords: Prediction regions System identification Nonlinear systems Uncertainty abstract This paper presents a new methodology to obtain prediction regions of the output of a dynamical system. The proposed approach uses stored past outputs of the system and it is entirely data-based. Only two hyperparameters are necessary to apply the proposed methodology. These scalars are chosen so that the size of the obtained regions is minimized while fulfilling the desired empirical probability in a validation set. In this paper, methods to optimally estimate both hyperparameters are provided. The provided prediction regions are convex and checking if a given point belongs to a computed prediction region amounts to solving a convex optimization problem. Also, approximation methods to build ellipsoidal prediction regions are provided. These approximations are useful when explicit descriptions of the regions are necessary. Finally, some numerical examples and comparisons for the case of a non-linear uncertain kite system are provided to prove the effectiveness of the proposed methodology. ©2023 The Author(s). Published by Elsevier Ltd on behalf of ISA. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 1. Introduction Dynamical models, either linear or nonlinear, suffer from a certain degree of modeling errors and uncertainty that results in inaccurate predictions of the outputs. One way to consider these uncertainties is by using prediction regions in which the future outputs are guaranteed to be contained with a pre-specified probability [1,2]. There are many applications where this kind of strategy is useful. For example, in stochastic programming and chance constrained problems [3,4], some of the parameters are random variables and, thus, it is necessary to characterize such uncertainty. Also, stochastic model predictive control (MPC) and chance constrained MPC [5] require some sort of uncertainty or disturbance characterization within the model of the system in order to solve the control problem. For the case of single-output systems, interval predictions can be used. This is a particular and easier case of the aforementioned prediction regions. These intervals can be constructed in two different ways, one consisting on following a robust scheme, that is, the prediction will be included in the interval with a probability equal to one. For this uni-variate setting, there are state-of-the- art techniques that work reasonably well. For example, if the uncertainty is bounded, set membership methods can be used to obtain robust intervals [6,7]. However, robust approaches can be overly conservative and thus a second approach has gained the ∗Correspondence to: Avenida Camino de los descubrimientos s/n, 41092, Sevilla, Spain. E-mail address: [email protected] (A.D. Carnerero). attention of researchers. This approach is the use of stochastic methods. In a stochastic setting, the prediction belongs to the aforementioned interval with a certain probability 1 −τ. In this context, one could rely on parametric regression techniques, obtaining a quantile regressor (QR) [8,9]. The main limitation of these methods is that a large number of training samples is necessary when τis close to the extremes of (0,1). Similarly, as an extension of QR techniques, one could rely on Interval Predictor Model (IPM) approaches [10]. The aforementioned methods can deal with prediction intervals, i.e., single output systems. However, computing prediction regions for multi-variate systems becomes a harder task. The simplest way would be to obtain intervals considering each variable independently and then construct rectangular prediction regions (see Section 2.2.3 in [11]). The main advantage of this method lies in its simplicity. However, the desired probability may not be attained or the size of the regions may be too large. Bootstrap methods are used in [12,13] to construct regions that contain a path of a random variable with at least a certain probability. This means that the obtained regions are actually intervals for ap−step prediction. In the case of multi-variable dynamical systems, the literature is rather scarce. For example, in [14], prediction regions for a simple multivariate linear regression model are obtained. On the other hand, [15] proposes a method to calculate prediction regions of a system by computing the Jacobian of the Partial Least-Squares Regression (PLS) parameters. Thus, by means of this local linearization, an ellipsoid-shaped region can be obtained. Also, [16] manages to obtain ellipsoid regions for dynamical systems by means of an Inverse Regression https://doi.org/10.1016/j.isatra.2023.03.048 0019-0578/©2023 The Author(s). Published by Elsevier Ltd on behalf of ISA. This is an open access article under the CC BY license (http://creativecommons.org/ licenses/by/4.0/). A.D. Carnerero, D.R. Ramirez, S. Lucia et al. ISA Transactions 139 (2023) 49–59 (IR) scheme. This IR scheme is more efficient and reliable than the classic regression approaches, when high dimensional data is available. In [17], a data-driven framework to generate and evaluate ellipsoidal prediction regions to characterize the uncertainty of a time-series is proposed. This methodology is applied to the electricity prices as a path forecasting problem. In the field of machine learning, Conformal Prediction techniques [18] are used to obtain prediction regions for any method producing a central prediction of the outputs. Similarly, it is proposed in [19] a framework to obtain a stochastic model based on a deterministic model of a robotic dynamic system. On the other hand, assuming that any finite set of samples follows a multivariate normal distribution, Gaussian processes (GPs) [20] can also be used. The main limitation of this approach is that it relies heavily on the knowledge of the first two moments of the underlying multivariate probability distribution. In this paper, we propose a new methodology to compute prediction regions of a multivariate dynamical system by means of dissimilarity functions. These functions have been used in [21,22] to obtain uni-variate interval predictors. In those papers, an empirical probability distribution is obtained and then the intervals are computed from the quantiles of such empirical distribution. Here, the dissimilarity functions are used in a different setting to obtain prediction regions for the multidimensional case. The proposed region is defined implicitly through a convex optimization problem, skipping the computation of the empirical distribution, that would be intractable for the multivariate case. Thus, the procedure is much lighter in terms of computational burden than the one presented in [21]. Also, the hyperparameters in which the methodology relies are computed to minimize the size of the regions while achieving the desired empirical probability. Moreover, the proposed approach does not make any assumption on the data or the shape of the region, being able to represent more tightly different types of regions. Also, in order to tackle the case where explicit regions are necessary, a method to obtain ellipsoid approximations of the regions is provided. The rest of the paper is organized as follows. Section 2presents the dissimilarity function that will be used throughout the manuscript whereas its role in forecasting the output of a nonlinear system is presented in Section 3. Section 4defines the predictions regions and proposes a method to tune the hyperparameters. Also, in Section 5an approximation method based on a quadratic upper bound of the optimization problem is proposed. By means of this approximation, an ellipsoidal region can be obtained. In Section 6, two numerical examples of nonlinear dynamical systems are proposed in order to compare the performance of the region predictor with other baselines. Finally, Section 7presents the conclusions. 2. Dissimilarity functions Dissimilarity functions [21] are used to measure the similarity of an element to the elements of a given set. They will play a key role in the remaining sections of the paper, thus a brief summary of the key concepts is given in the following. Consider a set D=[d1d2. . . dN]∈Rn×N,where di∈Rn. We are interested in determining how similar a given vector d∈Rnis with respect to the other vectors of the set D. Thus, a dissimilarity function J(·):Rn×D→ [0,∞], measures the dissimilarity between a given point dand the data set D. Large values of J(d,D) imply a high dissimilarity whereas small values correspond to a high similarity. A first candidate could be the minimum distance of the point dto each member in the data set D. That is J(d,D)=min i=1,...,N∥d−di∥,(1) Fig. 1. Elements of set D⊂R2(in red) and three points in R2(in blue). Note that d1 iand d2 idenote the first and second component of each direspectively. where ∥·∥is a given norm. Another possibility would be to consider the mean value of the distances of dto each member of D. However, these dissimilarity functions, although simple, are not invariant with respect to affine transformations and thus they depend on the choice of the coordinate system. A better dissimilarity function would be Jγ(d,D)=min λ1,...,λN (1 −γ) N ∑ i=1 wiλ2 i+γ N ∑ i=1 |λi|(2a) s.t. d= N ∑ i=1 λidi(2b) 1= N ∑ i=1 λi,(2c) where γ∈ [0,1) is a tuning parameter and wi>0 constant weights. This dissimilarity function returns a small value when the vector is very similar to the elements of Dand it returns a large value when the vector is very dissimilar. To illustrate this concept, consider Fig. 1 in which a set D(red markers) and three element examples (blue markers) denoted as d1,d2and d3are plotted. The values of the dissimilarity function with γ=0 and wi=1,∀ifor each regressor are: J0(d1,D)=0.0021,J0(d2,D)=0.0085,J0(d3,D)=0.0226. Thus, d1attains the lowest dissimilarity, i.e. it is more similar with respect to the set because the concentration of points around d1 is larger than around d2and d3. On the other hand, d3attains the highest dissimilarity due to the fact that d3is actually outside the region in which all the elements of Dlie. This optimization problem is strictly convex and thus it has a unique solution. In order to guarantee that this solution exists, we assume that [d1. . . dN 1. . . 1]is full row rank. To justify that problem (2) is strictly convex, consider that it can be rewritten as Jγ(d,D)=min λ(1 −γ)λ⊤Hλ+γ∥λ∥1(3a) s.t. [D 1]λ=[d 1],(3b) where ∥λ∥1=∑N i=1|λi|and H>0 because the weighting factors wiare strictly positive. It is easy to see that the quadratic term 50 A.D. Carnerero, D.R. Ramirez, S. Lucia et al. ISA Transactions 139 (2023) 49–59 Fig. 2. Relation between past inputs, outputs, the operator h(·) and the regressor zk. (1 −γ)λ⊤Hλis strictly convex because H>0 and (1 −γ)>0. Also, the 1-norm term ∥λ∥1is convex as well because γ∈ [0,1). Thus, the sum of the quadratic term and the absolute value term is a strictly convex function. Therefore, the optimization problem is strictly convex. Remark 1. This dissimilarity function is invariant to affine transformations, being more convenient than the ones based only on the distance of the point dwith respect to the data set D because they are dependent with respect to the reference system (see [21,23]). 3. Forecasting nonlinear systems using dissimilarity functions Consider a discrete time nonlinear system xk+1=f(xk,uk) yk=g(xk, vk),(4) where kis the time instant, xk∈Rnxare the states, uk∈Rnu are the inputs, yk∈Rnyare the outputs, vk∈Rnvaccounts for measurement noise, f(·) and g(·) are unknown nonlinear functions such that f(·):Rnx×nu→Rnxand g(·):Rnx×nv→Rny. We assume that only past outputs and inputs are accessible. In this section, we aim to obtain the prediction of the output value denoted as ˆ yk. This prediction will be computed as a function of a regressor denoted as zkand a given data set of past regressors and their corresponding output values. As we assumed that only past outputs and inputs are measurable, the regressor zkis obtained by applying a certain operator h(·) to past outputs and inputs, with h(·):Rnunpu+nynpy →Rnzwhere npu and npy are the number of past terms of the input and the output respectively, and nzis the dimension of the regressor (see Fig. 2). It is also assumed that this operator is known. The operator, included for generality, could be the identity, which would provide the raw past values of the inputs and the outputs without doing any mathematical operation, that is, zk= [uk−1⊤,...,uk−npu ⊤,yk−1⊤,...,yk−npy ⊤]⊤.However, more complex operators can be used, such as kernel or neural networks based operators (like in the numerical example of Section 6). Denoting by ¯ zi∈Rnzand ¯ yi∈Rnythe past samples of the regressor and the output respectively, we assume that the data sets are stored in the form of the matrices Z=[¯ z1¯ z2. . . ¯ zN],Y=[¯ y1¯ y2. . . ¯ yN], where Nis the number of samples. These data sets will be used to obtain the prediction of the system output for a given regressor z. Here, it is proposed to compute such prediction using a weighted sum of the elements of Y, that is, ˆ y(z)= N ∑ i=1 λ∗ i¯ yi.(5) We will discuss in the following how to compute the weights λ∗ i. To obtain the desired prediction of the output using the proposed structure, we consider the dissimilarity function given in the previous section written for the forecasting problem as Jγ([z y],[Z Y])=min λ1,...,λN (1 −γ) N ∑ i=1 wiλ2 i+γ N ∑ i=1 |λi|(6a) s.t. z= N ∑ i=1 λi¯ zi(6b) y= N ∑ i=1 λi¯ yi(6c) 1= N ∑ i=1 λi.(6d) We will compute ˆ y(z), the expected value for the output y given the value of the regressor z, as the one that gives the minimum value of the dissimilarity function, that is, ˆ y(z)=arg min yJγ([z y],[Z Y]).(7) The following lemma states how to find ˆ y(z) by solving a simpler optimization problem which in turn will return the optimal values of the weights necessary to compute ˆ y(z) using the structure of (5). Lemma 1 (Computation of the Prediction).The prediction ˆ y(z)can be obtained as: ˆ y(z)=Yλ∗= N ∑ i=1 λ∗ i¯ yi,(8) where λ∗is computed as λ∗=arg min λ1,...,λN (1 −γ) N ∑ i=1 wiλ2 i+γ N ∑ i=1 |λi|(9a) s.t. z= N ∑ i=1 λi¯ zi(9b) 51 A.D. Carnerero, D.R. Ramirez, S. Lucia et al. ISA Transactions 139 (2023) 49–59 1= N ∑ i=1 λi.(9c) Proof. Let us consider Jγ(z,Z), that is, Jγ(z,Z)=min λ1,...,λN (1 −γ) N ∑ i=1 wiλ2 i+γ N ∑ i=1 |λi|(10a) s.t. z= N ∑ i=1 λi¯ zi(10b) 1= N ∑ i=1 λi.(10c) As this optimization problem is a less constrained version of that in (6), the optimal value of the cost of the more constrained problem is always greater or equal to the optimal value of the cost of the less constrained one, that is, Jγ(z,Z)≤Jγ([z y],[Z Y]).(11) Now, denote the vector of weights λthat minimizes the optimization problem in (10) as λ∗. Then, it is easy to see that Jγ(z,Z)=Jγ([ z Yλ∗],[Z Y]). Because of this and (11), it is not possible to find a value of y improving the value of the dissimilarity function. Thus, ˆ y(z)= Yλ∗=∑N i=1λ∗ i¯ yi, computed using the optimal weights from (10), is the one that minimizes the dissimilarity function for a fixed z in Eq. (7).■ An important property is the uniqueness of the prediction. The following remark states that ˆ y(z) is unique. Remark 2 (Uniqueness Of ˆ y(z)).From the results stated in the previous section problem (6) is strictly convex and because [ZTYT1T]Tis assumed to be full rank, its solution exists. Therefore, for a certain z,Zand Y, the dissimilarity function (6) has a unique minimum at ˆ y(z)=Yλ∗, where λ∗corresponds to the optimal solution of problem (10), and therefore the prediction ˆ y(z) is unique. The numerical computation of problem (9) can be addressed by means of a dual formulation. For this particular optimization problem, the number of dual decision variables corresponds to the number of equality constraints, i.e. (nz+1) which is obviously much smaller than the number of data points (N). Also, note that once the dual variables are fixed, the primal variables can be obtained by solving Nunidimensional problems which have explicit solutions. In this paper, the numerical examples have been computed by using an accelerated gradient method in the dual variables [24,25]. For the univariate case (ny=1) and for the particular choice γ=0, it was proved in [21] that this forecasting corresponds to the one obtained by means of least-squares regression. As we explain in what follows, strictly positive values of γencourage the components of λto be positive (which means that the central estimation is often obtained from an interpolation of points). When every λi≥0,∀i=1,...,N,∥λ∥1=1. However, when some of the components of λare negative, we infer that ∥λ∥1>1. In other words, the cost term ∥λ∥1becomes larger when extrapolating points (i.e. using negative values of λi). This means that convex combinations of λiare encouraged and thus interpolation is preferred, which may improve the predictions. The choice of γwill be further discussed in the context of the proposed methodology in Section 4.4. Remark 3. The weighting factors wioffer a way to add flexibility to the forecasting scheme. They could be used to consider locality in the data, by making wi=max{∥z−zi∥2 2, ϵ}, where ϵis a small positive constant. In this way, local estimation can be implemented (see [21]). 4. Region estimation Besides obtaining a forecasting for the next output, the aforementioned dissimilarity function allows us to compute regions that will contain the output with a pre-specified probability. This is an extension of the method proposed in [21] where only one-dimensional outputs were considered and, thus, interval predictions were obtained. In the present work, we address the multidimensional case, providing prediction regions for the multidimensional output. The scheme proposed in [21] is based on a numerical integration in a unidimensional space. The generalization of [21] to the multidimensional setting is not trivial because of the well-known complexity of high-dimensional numerical integration. Therefore, that technique is not useful in practice for multivariable systems. In what follows, we show how to obtain prediction regions in multidimensional spaces without resorting to numerical integration, providing a fundamental advantage with respect to the results presented in [21]. In order to characterize such a region, two issues must be tackled: the choice of the region center and the computation of the region itself. 4.1. Choosing the center In the previous section, we denoted ˆ y(z) as the value of y that minimizes the dissimilarity function given z(7). This optimal value can be obtained as ˆ y(z)=Yλ∗where λ∗corresponds to the optimal solution of the strictly convex problem (10). As the prediction ˆ y(z) is unique (see Remark 2), it will be used to define the center of the region that will contain the real value of the predicted output with a given probability. From now on, in order to make the manuscript more readable, we simplify the notation by removing Zand Yfrom the dissimilarity function, since they are assumed to be fixed. That is, we use the notation Jγ([z y],[Z Y])=Jγ(z,y). Moreover, the minimal value of the dissimilarity function given z is denoted as J∗ γ(z). That is, J∗ γ(z)=Jγ(z,ˆ y(z)) =min yJγ(z,y). 4.2. Computing the regions The dissimilarity function can be seen as a sort of surrogate of the probability distribution of ygiven z[21]. Thus, the prediction region will be defined in relation to the value of J∗ γ(z). We consider that this probability distribution has a maximum at ˆ y(z) and decreases as the dissimilarity function increases. For that reason, we consider prediction regions that are defined as those points for which the dissimilarity function does not exceed more than a given factor αthe one corresponding to the central prediction ˆ y(z). This is formally stated in the following definition. Definition 1 (Prediction Region).For a given z,γ, data sets Z,Y and a tunable parameter α > 1, we define a prediction region as the set ∆(z)={y:Jγ(z,y)≤αJ∗ γ(z)},(12) that is, the points ythat obtain a dissimilarity less or equal to αJ∗ γ(z). 52 A.D. Carnerero, D.R. Ramirez, S. Lucia et al. ISA Transactions 139 (2023) 49–59 Table 1 Numerical results of the clarifying example. Area Empirical probability γ=0 1.1423 0.9600 γ=0.2 1.0077 0.9580 γ=0.4 0.9834 0.9640 γ=0.8 0.9561 0.9520 We note that the provided regions are defined in an implicit way (we address explicit ellipsoidal regions in the next section). Since the dissimilarity function is convex in y, the obtained implicit regions are convex and can be used in different settings, e.g. in chance-constrained optimization. For example, in many cases, it is not necessary to compute a prediction region, but to verify if a certain point ¯ ybelongs to it. This is an affordable task as it suffices to compute the dissimilarity Jγ(z,¯ y) and check that (12) holds. 4.3. Clarifying example: multivariate uniform distribution To illustrate the role of γ, consider the easier case of a multivariate uniform distribution. As we are not considering a dynamic system, conditioned distributions are not taken into account. A number of samples N=500 are extracted from a uniform distribution in R2and gathered into a data set Y. This uniform distribution is generated considering that both dimensions are independent one with respect the other, with values ranging from 0 to 1. Then, the obtained distribution is rotated with an angle of π 6rad in order to misalign the previously computed distribution. We will approximate the support of the generated distribution by means of prediction regions of the form ∆={y:Jγ(y,Y)≤αJ∗ γ}. Here, we consider a finite family of four possible values of γ, denoted as Γ= {0,0.2,0.4,0.8}. As we stated in the previous sections, γ=0 corresponds to the least-squares solution [21] and thus an ellipsoid will be obtained. However, more appropriate shapes can be obtained with different values of γ. Also, note that the empirical expectation turns out to be the center of the region. In this example, we consider τ=0.05, that is, we would like the probability of the samples falling within the region to be 95%. The results are shown in Fig. 3. The blue crosses correspond to the points falling inside the region whereas the red crosses correspond to the points falling outside. It is easy to see that γ=0.8 is the element of Γthat captures best the shape of the proposed probability distribution, leading to a tighter region. Also, it can be seen in Table 1 that, for this example, increasing the value of γprovides a smaller area (tighter approximation) while still fulfilling the desired empirical probability. Even though the difference might seem small at first, we note that this difference may increase when considering regions of greater dimensions. 4.4. Tuning the hyperparameters αand γ Regarding which value of γshould be used, consider a validation set V= {(˜ z1,˜ y1),(˜ z2,˜ y2),...,(˜ zNv,˜ yNv)}, composed of pairs (˜ zj,˜ yj) gathered from the system. The problem of determining the best value of γ∈ [0,1) can be addressed in different ways. Suppose that, in order to reduce the complexity of the problem, γis constrained to belong to a set Γ⊆ [0,1) of finite cardinality (e.g. a set of equidistant points in the interval [0,1)), i.e., a finite family of γ. Then, a reasonable choice would be the value that minimizes the error of the central predictions {ˆ y(˜ z)}with respect to the real outputs {˜ y}in the validation set V, that is: γ∗=arg min γ∈Γ 1 Nv Nv ∑ j=1 ˆ y(˜ zj)−˜ yj  2 2. Note that minimizing the errors of the predictions is related to minimizing the size of the regions. This is due to the fact that reducing the gap between ˆ y(˜ zj) and ˜ yjleads to a smaller value of αneeded to fulfill the desired probability 1 −τ. Given γ, smaller values of αmake the regions smaller. This can be easily seen by looking at Definition 1, where the equation of the prediction region is presented. There, we note that αappears multiplying the optimum value J∗ γ(z), i.e. the optimum value of the dissimilarity function given z. Therefore, αis a parameter that implicitly defines the size of the region. That is, for the case when α=1, we have that the region is composed of just the point ˆ y(˜ zj) whereas for larger values of α, the number of points included in the region increases. The procedure to obtain αpresented in the following aims to obtain the smallest possible region that guarantees that a point ˜ yjtaken from the validation set Vis contained in the computed region ∆(˜ zj) with a pre-specified probability 1 −τ. In order to circumvent the difficulty of generating i.i.d. samples {(z,y)}, we consider that the discrete probability described in the validation set Vis assumed to represent the true probability distribution. Then, αis chosen so that it fulfills ProbV(y∈∆(z)) ≥1−τ, that is, guaranteeing that the fraction of points yin the validation set that fall out of the corresponding prediction region does not exceed τ. Denote ˜αj=Jγ(˜ zj,˜ yj) J∗ γ(˜ zj),j=1,...,Nv. Then, the value of αsatisfying the probabilistic specification in the set Vis the one corresponding to the rτ= ⌈τNv⌉largest value of { ˜αj}Nv j=1. This choice makes the number of points that fall out of the prediction region equal to rτ−1 which implies that the fraction of such points out of the region is equal to rτ−1 NV= ⌈τNv⌉−1 Nv≤τNv Nv=τ. Thus, the number of points out of the region do not exceed τas desired. On the other hand, if i.i.d. samples are available, one could resort to the concept of probabilistic scaling to choose the value of rτand Nv(see [26–28]). In that case, two levels of probability results are obtained that apply to any pair (¯ z,¯ y) generated according to the underlying distribution of (z,y). 4.5. Comparison with the approach presented in [21] In this subsection, we show that the proposed approach is superior to the one presented in [21] in terms of computational complexity. For this purpose, we consider the example presented in [21] corresponding to the Lorenz attractor. This is a chaotic system following the set of differential equations: do dt=σ(p−o) dp dt=o(ρ−q)−p dq dt=op −βq, where σ,ρand βare real scalar parameters. In this example, these parameters take the values σ=10, ρ=28 and β=8/3. Considering a sample time of Ts=0.1 s, we tackle the problem of obtaining 1−step ahead interval predictions of the output of the system, that is yk=ok. The regressor zkis composed of the 53 A.D. Carnerero, D.R. Ramirez, S. Lucia et al. ISA Transactions 139 (2023) 49–59 Fig. 3. Prediction regions for different values of γ. Fig. 4. Computed intervals for the proposed approach (τ=0.2). two previous values of the outputs, i.e. zk= [yk−1,yk−2]⊤. The size of the data set, validation set and test set is 300, 500 and 500 respectively. The considered finite family of γ, that is Γ, is composed of the values [0, 0.1, 0.2]. Also, we consider the same grid on the probability distribution as in [21]. The results are shown in Fig. 4 and Table 2. There, it can be seen that the new approach attains a much smaller computational time. Also, the obtained interval widths and empirical probabilities are better (that is, the interval is smaller and the empirical probability higher). Note that the computational time corresponds to the total amount of required time to calculate the hyperparameters and compute the intervals for all the points in the test set. Table 2 Numerical results of the example in [21]. Interval width Empirical probability Time (s) Approach in [21] 2.0352 0.8900 296.47 Proposed 1.9773 0.8980 6.9201 5. Ellipsoidal regions Due to the term γ∑N i=1|λi|appearing in the definition of the dissimilarity function, the prediction regions cannot be computed explicitly. Instead, in order to check if a point ybelongs to a prediction region, a convex optimization problem must be solved. We provide in this section a quadratic upper bound to the dissimilarity function that leads to ellipsoidal regions that can be computed in an explicit way. The idea is to upper bound the absolute values {|λi|}N i=1with scalar quadratic functions. This transforms the functional into a quadratic one and thus the optimization problem can be solved explicitly since the minimization of a convex quadratic function subject to linear equality constraints has an explicit solution. In order to upper bound the absolute values, we use the following lemma. Lemma 2 (Quadratic Upper Bound Of |λi|).Given a scalar λc iand ν > 0, |λi| ≤ |λc i| + ν 2((|λi| |λc i| + ν)2 +1),∀λi∈R.(13) Proof. See Appendix B in [29]. ■ We choose λc ias the λ∗ iused to obtain the central prediction ˆ y(z)(8). We note that the upper bound provided by Lemma 2 is tight when λi=λc i=λ∗ iand ν→0. Thus, we choose νas a 54 A.D. Carnerero, D.R. Ramirez, S. Lucia et al. ISA Transactions 139 (2023) 49–59 small constant. Applying Eq. (13) to the optimization problem in (6), we obtain Jγ(z,y)≤Qγ(z,y), where Qγ(z,y)=min λ1,...,λN 1 2λTHγλ+cγ(14a) s.t. [Z 1]λ=[z 1](14b) Yλ=y,(14c) Hγis a diagonal matrix with entries Hγ|i,i=2(1 −γ)wi+γ |λc i| + ν,i=1,...,N and cγ=γ∑N i=1 |λc i|+ν 2. The following lemma characterizes the new dissimilarity function Qγ(·). Lemma 3. Assume that Hγ>0and that the matrix [ZTYT1T]Tis full rank. Denote Qγ(z,y)=min λ1,...,λN 1 2λ⊤Hγλ+cγ s.t. [Z Y 1]λ=[z y 1], y∗=arg min yQγ(z,y), Q∗ γ(z)=Qγ(z,y∗). Then, y∗=Γ1,2⊤Γ−1 1,1[z 1],Q∗ γ(z)=1 2[z 1]⊤Γ−1 1,1[z 1]+cγ Qγ(z,y)=Q∗ γ(z)+1 2(y−y∗)⊤Φ2,2(y−y∗). where Γ1,1=[Z 1]H−1 γ[Z 1]⊤,Γ1,2=[Z 1]H−1 γY⊤, Γ2,2=YH−1 γY⊤,Φ2,2=(Γ2,2−Γ1,2⊤Γ−1 1,1Γ1,2)−1. Proof. The proof is shown in Appendix.■ We now define, by means of the new dissimilarity function Qγ(z,y), the following ellipsoidal regions Eγ ,α(z)={y:Qγ(z,y)≤αQ∗ γ(z)} ={y:(y−y∗)⊤Φ2,2(y−y∗)≤2(α−1)Q∗ γ(z)}, where, the values for y∗,Φ2,2and Q∗ γ(z) are given in Lemma 3. We note that, given γ≥0, the optimal value for αcould be obtained by means of a validation set Vas in the case of implicit regions. In this case, the scalars ˜αjare defined as ˜αj=Qγ(˜ zj,˜ yj) Q∗ γ(˜ zj),j=1,...,Nv. The optimal value for αwould be the rτlargest one, where rτ= ⌈τNv⌉. We also remark that, in this case, one could obtain γin such a way that a measure of the size of the obtained ellipsoids is minimized. 6. Numerical example In this section, we propose as an example the application of the proposed methodology to compute prediction regions for the predictions of the outputs of two different multivariable systems. 6.1. Towing kite In this section, we propose as an example the application of the proposed methodology to compute prediction regions for the predictions of the outputs of a multivariable system. As the system to be forecasted, we consider the non-linear model of a towing kite presented in [30]. It has three states, which corresponds to the angles of the kite, θ,φand ϕ, a control input uand two uncertain parameters. The equations governing the system are the following ˙ θ=va LT(cos ϕ−tan θ E), ˙ φ= − va LTsin θsin ϕ, ˙ϕ=va LT u+˙ φcos θ, where va=v0Ecos θ, E=E0−0.028u. From this system of equations, the parameter LT(length of the tether) is considered to be known whereas v0(wind speed) and E0 (base glide ratio) are considered to be uncertain parameters. Note that these uncertainties affect directly va, which can be described as the effect of the wind and the glide ratio E, and thus it also affects the states θ,φand ϕ. Furthermore, only the states θand φare considered as outputs with an additive correlated gaussian noise so that N(µ, Σ):µ=[0 0],Σ=[0.01 −0.0085 −0.0085 0.01 ]. We assume that the control input uis computed following a certain stabilizing control policy κ(θ, φ, ϕ) whose exact nature is irrelevant for this paper (see [31] for more details on a possible control law). In fact, we treat the kite as an autonomous closed loop system. Thus, it is assumed that we have some stable trajectories of the closed-loop system for different values of the base glide ratio E0and the wind speed v0that are used to build the data sets Z,Ywith a sample time of Ts=0.15 s. We consider two different data sets in this example, one with N=250 samples and another one with N=500 samples. Also, a validation set V composed of Nv=500 samples to compute the optimal values of γand αis considered. On the other hand, we have a test set Tof NT=1000 samples to compare the obtained results with other baseline strategies. We consider a finite family of γdenoted as Γ. This Γis composed of the values Γ= {0,0.2,0.4,0.6,0.8}. In this example, the operator h(·) is chosen as a Neural Network (NN). That is, the NN is working as an encoder network. This means that the feature space of the whole data sets Zand Yare transformed by means of this given NN. In this case, the NN is comprised of 2 hidden layers of 10 neurons each one. Furthermore, the inputs of the NN are the last 5 outputs of the system. The outputs of the last layer are used as the regressor z. Note that this regressor zis the same for the proposed approach and the baselines and thus the following numerical comparisons are fair. Also, it is assumed that this NN is given. This means that the problems related to training appropriately the NN, finding the optimal number of layers, etc. are outside the scope of this paper. 55 A.D. Carnerero, D.R. Ramirez, S. Lucia et al. ISA Transactions 139 (2023) 49–59 Fig. 5. Obtaining the data sets Zand Yfor the kite example. Note that yk−jcorresponds to (θk−j, φk−j). Fig. 6. Projected regions for a representative trajectory (τ=0.2). Fig. 5 shows the flow of the data used to build the data sets Z and Yin this example. There, we note that the previous values of the output go into the neural network, whose output is the regressor zk. On the other hand, the previous state along with the control input given by the control law κ(θk−1, φk−1, ϕk−1) go into the system, producing the next value of the output yk. Then, the pair (zk,yk) is stored in the data sets Zand Y. Three baseline strategies are considered for the sake of comparison. First, a quantile regression (QR) approach. This technique is well established for the single output case. In order to obtain regions, a probabilistic interval of probability τis obtained for each output independently and then combined to form a boxshaped region. However, combining two intervals of probability τdoes not turn out to be a prediction region of probability τ. Actually, a lower probability will be attained. To tackle this problem, a back-off parameter is calculated so that the computed boxes fulfill the empirical probability in the validation set V. The second one is based on GPs. Here, the matlab function ‘‘fitgrp’’ is used to obtain a model for each output independently. This function automatically finds an appropriate Kernel function for the input–output data and optimizes the values of the hyperparameters. Then, using ‘‘predict’’, it is possible to obtain intervals of a specified probability τ. Same as before, it is necessary to compute a back-off parameter in the validation set Vto achieve the desired empirical probability. Finally, an implementation of the prediction regions based on Inverse Regression (IR) of [16] is considered. Fig. 6 shows the obtained predictions projected as intervals using the proposed approach for the case of a representative trajectory of the test set employing τ=0.2 and N=250. It can be seen how the noisy trajectory (in red) is enclosed by the computed region most of the times, as it was expected. The results for τ=0.1 and τ=0.2 are shown in Table 3. It can be seen how the implicit regions of the proposed approach and the approximation ellipsoids outperform the baselines considered for both values of τand different data set sizes. That is, the proposed approach obtains regions of considerably smaller area while still fulfilling the specified probability given by τ(i.e. with an empirical probability (E.P.) no smaller than 1−τ). This means that the proposed approach provides regions of smaller uncertainty in comparison to the baselines considered. It is clear that obtaining smaller regions will allow us to take better decisions if we implement the proposed methodology in an optimization setting. The average time required for N=250 to check if a certain point belongs to the region is 5.813 ms for the proposed approach (we allude here to the true region not the approximation). On the other hand, for the case of GPs, QR and IR, the values are 0.011 ms, 0.001 ms and 0.023 ms respectively. 56 A.D. Carnerero, D.R. Ramirez, S. Lucia et al. ISA Transactions 139 (2023) 49–59 Table 3 Area of the regions and empirical probabilities obtained for the proposed methodology, the approximation with ellipsoidal regions, Gaussian processes (GPs), quantile regression (QR) and inverse regression (IR). Proposed Approximation GPs QR IR N Area E.P. Area E.P. Area E.P. Area E.P. Area E.P. τ=0.1250 0.1695 0.9200 0.1885 0.9170 0.2235 0.9250 0.5839 0.9120 0.4331 0.9031 500 0.1551 0.9260 0.1656 0.9170 0.2092 0.9240 0.2645 0.8960 0.4101 0.9043 τ=0.2250 0.1119 0.8070 0.1176 0.8070 0.1518 0.8450 0.2341 0.8550 0.2466 0.7446 500 0.1053 0.8150 0.1142 0.8260 0.1403 0.8310 0.1590 0.8080 0.2335 0.7454 Fig. 7. Projected regions for a representative trajectory (τ=0.2). 6.2. Rössler attractor As a different example, consider the system described by the following differential equations do dt= −p−l dp dt=o+ap dl dt=b+l(o−c), which is known as the Rössler attractor. In this example, the parameters take the values a=0.2, b=0.2, c=5.7 and the sample time is 0.1 s. The aforementioned system along with the proposed parameters show a chaotic behavior. Two different sizes (150 and 250 samples) for the data set (Zand Y) are available. For the validation set Vand the test set T, we consider a total of 1000 samples in each one. The regressor is comprised of the past output of the system, i.e. zk= [ok−1,pk−1]⊤. Note that this corresponds to choosing h(·) as the identity operator. Also, the output of the system is affected by some gaussian noise so that N(µ, Σ):µ=[0 0],Σ=[0.05 0.025 0.025 0.05 ]. The baselines used for comparison will be the same as in the previous example: GPs, QR and IR. However, in this example, the IR method fails to obtain the desired regions (the obtained regions are several orders of magnitude larger), leading to unacceptable results. For this reason, we removed IR from the table of results. The results are shown in Fig. 7, and Table 4.Fig. 7 shows the real evolution of the output and the prediction region obtained with the proposed approach and projected as in the previous example. Table 4 summarizes the numerical results of the experiment. It can be seen that the proposed strategy achieves the smallest areas. On the other hand, the ellipsoidal approximation obtains similar areas with respect to GPs when considering the small data set (150 samples) and smaller areas when considering the larger data set (250 samples). 7. Conclusions This work presents a new data-based approach to obtain prediction regions to be used in nonlinear systems. The hyperparameters in which the methodology relies are computed so that the empirical probability is at least higher than the desired one in a validation set. Also, it is shown that checking if a point belongs to the computed region is as simple as solving a convex optimization problem. Finally, the proposed method has been used to generate prediction regions of a non-linear system, which have been compared favorably with the ones obtained by means of other baselines strategies based on Gaussian processes, quantile regression and inverse regression. As a future work, we consider the application of the proposed strategy in the context of developing novel stochastic or chance constrained MPC approaches, where the output of the prediction model is a whole region which heavily depend on uncertainty quantification and therefore could leverage the results of this paper. Another field of application would be failure detection, where the proposed technique could be used to obtain the regions of normal operation of the system. These regions could be then used to detect abnormal behaviors, which could be attributed to failures of the system. 57