Probabilistic data-driven methods for forecasting, identification and control

Escuela T´ecnica Superior de Ingenier´ıa Departamento de Ingenier´ıa de Sistemas y Autom´atica Doctoral Thesis Probabilistic data-driven methods for forecasting, identification and control Alfonso Daniel Carnerero Panduro Supervised by: Daniel Rodr´ıguez Ram´ırez Teodoro ´ Alamo Cantarero Seville, September 2022 Contents Acknowledgements i Abstract iii Notation, conventions and definitions v 1 Introduction 1 1.1 Motivation and objectives . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 The problem of system identification . . . . . . . . . . . . . 1 1.1.2 Quantifying the uncertainty . . . . . . . . . . . . . . . . . . 4 1.1.3 Model predictive control . . . . . . . . . . . . . . . . . . . . 7 1.1.4 Objectives of this dissertation . . . . . . . . . . . . . . . . . 11 1.2 Outline ................................. 11 1.3 Publications............................... 13 I Probabilistic forecasting 15 2 Forecasting using dissimilarity functions 17 2.1 Proposed dissimilarity function . . . . . . . . . . . . . . . . . . . . 18 2.1.1 Clarifying example . . . . . . . . . . . . . . . . . . . . . . . 22 2.2 Dissimilarity functions and regression . . . . . . . . . . . . . . . . 22 2.3 Application: forecasting stock prices using dissimilarity functions . 26 2.3.1 Results ............................. 27 2.4 Conclusions............................... 31 3 Probabilistic prediction regions 33 3.1 Univariatecase ............................. 33 3.1.1 Empirical probability density function . . . . . . . . . . . . 34 3.1.2 Clarifying example: uniform distribution . . . . . . . . . . . 38 3.1.3 Numerical example: Lorenz attractor . . . . . . . . . . . . . 38 3.1.4 Numerical example: Dow Jones industrial average index . . 40 3.2 Multivariatecase............................ 44 3.2.1 Implicit regions . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.2.2 Clarifying example: multivariate uniform distribution . . . 46 3.2.3 Ellipsoidal prediction regions . . . . . . . . . . . . . . . . . 46 3.2.4 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3 Conclusions............................... 52 II Kriging-based identification 53 4 State-space kriging for autonomous systems 55 4.1 Dynamickriging ............................ 55 4.2 Linear state-space kriging . . . . . . . . . . . . . . . . . . . . . . . 57 4.2.1 Initial condition . . . . . . . . . . . . . . . . . . . . . . . . 59 4.2.2 Local-data approach . . . . . . . . . . . . . . . . . . . . . . 59 4.3 Kernel-based state-space kriging . . . . . . . . . . . . . . . . . . . 60 4.3.1 Initial condition . . . . . . . . . . . . . . . . . . . . . . . . 63 4.4 Kalman filter for SSK . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.5 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.5.1 Sunspot number . . . . . . . . . . . . . . . . . . . . . . . . 65 4.5.2 R¨ossler attractor . . . . . . . . . . . . . . . . . . . . . . . . 66 4.6 Conclusions............................... 68 5 State-space kriging for non-autonomous systems 69 5.1 Non-autonomous linear SSK . . . . . . . . . . . . . . . . . . . . . . 69 5.1.1 Initial condition . . . . . . . . . . . . . . . . . . . . . . . . 71 5.2 Non-autonomous kernel-based SSK . . . . . . . . . . . . . . . . . . 71 5.2.1 Initial condition . . . . . . . . . . . . . . . . . . . . . . . . 72 5.3 ApplicationtoMPC .......................... 73 5.3.1 Nominal stability analysis . . . . . . . . . . . . . . . . . . . 75 5.3.2 Robust stability analysis . . . . . . . . . . . . . . . . . . . . 76 5.4 Examples ................................ 77 5.4.1 Continuously-stirred tank reactor . . . . . . . . . . . . . . . 78 5.4.2 Temperature control lab . . . . . . . . . . . . . . . . . . . . 79 5.5 Conclusions............................... 81 III Probabilistically-certified data center management 87 6 Bounds on the constraint violation level 89 6.1 Introduction............................... 89 6.2 Mainresults............................... 90 6.2.1 A first bound on constraint violation rate . . . . . . . . . . 91 6.2.2 A different bound . . . . . . . . . . . . . . . . . . . . . . . 92 6.3 ClarifyingExample........................... 94 6.4 Conclusions............................... 97 7 Energy-efficient management of data centers 99 7.1 Introduction............................... 99 7.2 Data center description . . . . . . . . . . . . . . . . . . . . . . . . 101 7.2.1 Tasksmodel...........................102 7.2.2 Servermodel ..........................102 7.2.3 Thermalmodel.........................104 7.2.4 Quality of service . . . . . . . . . . . . . . . . . . . . . . . . 106 7.3 Management approach . . . . . . . . . . . . . . . . . . . . . . . . . 106 7.4 Particle based solvers for complex optimization problems . . . . . . 108 7.4.1 Scenario-based approach . . . . . . . . . . . . . . . . . . . . 111 7.4.2 Parallel implementation . . . . . . . . . . . . . . . . . . . . 113 7.5 Bounds on the constraint violation rate . . . . . . . . . . . . . . . 114 7.6 Numericalresults............................ 114 7.6.1 QoS violation rate . . . . . . . . . . . . . . . . . . . . . . . 116 7.6.2 Thermal constraint violation rate . . . . . . . . . . . . . . . 117 7.6.3 Parallel computation improvement . . . . . . . . . . . . . . 118 7.6.4 Computation time analysis . . . . . . . . . . . . . . . . . . 119 7.7 Conclusions...............................120 8 Conclusions and future work 121 8.1 Contributions..............................121 8.2 Futurework...............................122 Bibliography 125 Acknowledgements Me voy a tomar la libertad de escribir esta secci´on en castellano debido a que algunas de las personas a las que va dirigida no tienen un conocimiento fluido de la lengua de Shakespeare y, por tanto, todo este esfuerzo podr´ıa resultar en vano. Adem´as que, despues de escribir numerosas p´aginas en riguroso ingl´es, no est´a de m´as volver a las propias ra´ıces por un tiempo. Han pasado aproximadamente 9 a˜nos desde que entr´e a estudiar en la escuela de ingenieros y, hasta cierto punto, parece como si fuera ayer (lo cual es probablemente el t´opico m´as usado de la historia). El c´omo llegu´e aqu´ı me parece, en cierto sentido, pura casualidad, pues durante mi m´as tierna infancia y adolescencia, mis intereres laborales y profesionales variaron con frecuencia. En cualquier caso, creo que debo gran parte de mi inter´es a la ingenier´ıa y al control a Antonio Nuevo y a Juanma Esca˜no. Si no hubiera frecuentado aquel taller de automatizaci´on en el Colegio Altair cuando me encontraba en Primaria, probablemente mi trayectoria hubiera sido distinta. Por otro lado, debo a mi hermano Jos´e Mar´ıa mi inter´es por la carrera investigadora, la cual ´el comenz´o mucho antes que yo en otra facultad de esta misma universidad. Por alg´un motivo, pens´e que se trataba de un trabajo fascinante con una calidad de vida bastante buena. En otro orden de cosas, mi padre result´o ser un gran apoyo durante mis primeros pasos en la escuela, debido a las experiencias similares que hab´ıa sufrido durante su juventud. En cuanto a mis hermanas, Ana y Guadalupe, les debo el no vestir como un pordiosero (lo cual no es poco). Bueno, algo m´as habr´a, pero eso es lo primero que se me ha venido a la mente. Para mi madre, s´olo se me ocurre decir que es ´unica en todos los sentidos. Tambi´en me gustar´ıa agradecer especialmente a Sergio Luc´ıa por su acogida en la Universidad T´ecnica de Dortmund durante mi estancia de doctorado. Fue una de esas experiencias positivas en las que uno descubre que el resto del mundo no es como Sevilla (quien tenga o´ıdos que oiga). Por supuesto, y como no pod´ıa ser de otra manera, tambi´en agradezco la labor que han tenido mis directores Dani y Teo en este trabajo y por la formaci´on recibida a lo largo de todo este tiempo. Aunque nuestros caminos se separen temporalmente, espero que vuelvan a juntarse en un futuro cercano. Tambi´en, me gustar´ıa agradecer especialmente a Dani Lim´on por la oportunidad que me brind´o en su momento. Gracias a eso, pude empezar a trabajar en este grupo de investigaci´on. A su vez, me gustar´ıa agradecer a los compa˜neros de doctorado que me ayudaron en mis comienzos, especialmente a Pablo y Pepe. De manera similar, tengo tambi´en un hueco particular para mis compa˜neros de caf´es y/o desayunos Joaqu´ın y Jos´e Antonio. Por otra parte, existen a su vez gran cantidad de personas que han contribuido i infinitesimalmente a convertirme en lo que soy a d´ıa de hoy y a las que, por tanto, estoy muy agradecido. No quiero realizar una enumeraci´on exhaustiva porque eso resulta sumamente peligroso (la probabilidad de olvidar a alguien es bastante alta y no es directamente proporcional a la importancia de la persona en cuesti´on). Por tanto, si est´as leyendo esta aburrida secci´on, puedes concluir sin temor a equivocarte que eres una de esas personas. Por ´ultimo, y no por ello menos importante, me gustar´ıa dedicar esta tesis a mis otros dos hermanos, Juanma y Alberto, los cuales no creo que tengan un inter´es especial por esta dedicatoria. En cualquier caso, ah´ı queda dicho. Alfonso Daniel Carnerero Panduro, Sevilla, septiembre de 2022. ii Abstract This dissertation presents contributions mainly in three different fields: system identification, probabilistic forecasting and stochastic control. Thanks to the concept of dissimilarity and by defining an appropriate dissimilarity function, it is shown that a family of predictors can be obtained. First, a predictor to compute nominal forecastings of a time-series or a dynamical system is presented. The effectiveness of the predictor is shown by means of a numerical example, where daily predictions of a stock index are computed. The obtained results turn out to be better than those obtained with popular machine learning techniques like Neural Networks. Similarly, the aforementioned dissimilarity function can be used to compute conditioned probability distributions. By means of the obtained distributions, interval predictions can be made by using the concept of quantiles. However, in order to do that, it is necessary to integrate the distribution for all the possible values of the output. As this numerical integration process is computationally expensive, an alternate method bypassing the computation of the probability distribution is also proposed. Not only is computationally cheaper but it also allows to compute prediction regions, which are the multivariate version of the interval predictions. Both methods present better results than other baseline approaches in a set of examples, including a stock forecasting example and the prediction of the Lorenz attractor. Furthermore, new methods to obtain models of nonlinear systems by means of input-output data are proposed. Two different model approaches are presented: a local data approach and a kernel-based approach. A kalman filter can be added to improve the quality of the predictions. It is shown that the forecasting performance of the proposed models is better than other machine learning methods in several examples, such as the forecasting of the sunspot number and the R¨ossler attractor. Also, as these models are suitable for Model Predictive Control (MPC), new MPC formulations are proposed. Thanks to the distinctive features of the proposed models, the nonlinear MPC problem can be posed as a simple quadratic programming problem. Finally, by means of a simulation example and a real experiment, it is shown that the controller performs adequately. On the other hand, in the field of stochastic control, several methods to bound the constraint violation rate of any controller under the presence of bounded or unbounded disturbances are presented. These can be used, for example, to tune some hyperparameters of the controller. Some simulation examples are proposed in order to show the functioning of the algorithms. One of these examples considers the management of a data center. Here, an energy-efficient MPC-inspired iii 4Chapter 1. Introduction Combinations of some of the aforementioned methods have been explored as well, i.e. neuro-fuzzy methods [37, 38] combining Takagi-Sugeno models and Neural Networks. Recently, another technique that is getting increasing attention is the Koopman operator [39, 40], where the nonlinear dynamics of the system are converted into linear dynamics in a infinite-dimensional state vector. Although the system becomes infinite-dimensional, it is possible to approximate them by a sufficiently large state-vector and apply methods from linear control theory. All of these methods assume that the dynamics of the system are completely unknown. However, it may exist situations where some partial knowledge of the system can be used alongside some black-box model. This leads us to the hybrid modeling framework [41, 42]. Here, the outputs are computed as the sum of two terms, one corresponding to the known dynamics and another one corresponding to the unknown dynamics. Note that this is very different to the grey-box modeling approach mentioned before. Figure 1.2: From first-principles models to black-box models. 1.1.2 Quantifying the uncertainty However, sometimes we are interested not only in obtaining nominal predictions of a certain system, but also a quantification of their uncertainty. That is, when the considered system presents noisy measurements and/or additive disturbances, it would be useful to characterize a region where the real output of the system may be with a specified probability instead of only computing the expected value. First, we focus on the univariate case, i.e. interval predictions. Given the regressor rk, the objective is to compute an interval I(rk)=[y− k, y+ k] such that we maximize the probability that ykbelongs to I(rk) while minimizing the interval width (y+ k− y− k). These two conflicting objectives can be reconciled if one minimizes the interval width subject to the constraint that I(rk) contains ykwith a pre-specified probability. Interval predictions play a relevant role in the control of uncertain systems. Zonotopes and DC Programming are used to obtain interval state estimators in [43] and [44] respectively. Interval observers for linear time-varying systems have been proposed in [45] and [46]. Fault detection methods based on zonotopic bounds 1.1. Motivation and objectives 5 can be found in [47]. In [48], set theoretic approaches are also used in the context of fault detection. Set membership methods [49, 50] can also be used to obtain interval predictions. A mixed Bayesian/set-membership approach is proposed in [51]. There exists different methods in the literature that address the problem of obtaining interval predictions. For example, if the vector of uncertainty is bounded and the considered system satisfies some Lipschitz assumptions, one can resort to bounded error methods [52] that guarantee that ykis always contained in I(rk). See, for example, [53] and [27]. Other bounded error strategies have been proposed in [54, 55, 56]. The statistical characterization of noise and disturbances can be used to enhance the performance of interval estimation methods. See [57, 58, 59] and references therein. Also, probabilistic validation methods can be used to assess the performance of the interval predictors [60, 61, 62, 63]. An important concept is that of quantiles [64, 65]. Denote Fyk(a|rk) as the cumulative distribution function of the associated output ykconditioned to the regressor rk, that is, Fyk(a|rk) = Prob{yk≤a|rk}, where ais a scalar. Given rk, we say that aτis the conditioned τ-quantile if Fyk(aτ|rk) = Prob{yk≤aτ|rk}=τ. The estimation of the conditioned quantiles is relevant in multiple applications (see [66] and [67]) and can be addressed using different methodologies. The most classical approach relies on the assumption that ykand rkare jointly normal. That is, the assumption that the (joint) probability density function of the (random) variables yand ris a multivariable normal probability density function. Under this assumption, the conditioned p.d.f. is a monovariable normal p.d.f. and the quantiles can be obtained in a simple and direct way [68]. Unfortunately, the methods based on normal distributions are very sensitive to the presence of outlier contamination. Moreover, in many long-tailed distributions, the normal assumption is not well suited to characterize confidence intervals and one has to resort to non-Gaussian distributions. In these cases, generalizations of the Chebyshev inequality can be used to obtain probabilistic bounds [69, 70]. The computation of the conditioned quantiles can be also addressed by means of parametric regression techniques [65], [66]. Assuming that there exists θfor which yk≈θ⊤rk, then parameter vector θcan be chosen as the one that minimizes a cost function of the error θ⊤rk−yk. If one chooses a cost function that penalizes in an asymmetric way positive and negative errors then a quantile regressor is obtained. Given the training pairs (¯yi,¯ri), i= 1, . . . , N and τ∈(0,1), the quantile regressor is defined in terms of the following optimization problem min θ N X i=1 (1 −τ) max{0, θ⊤¯ri−¯yi}+τmax{0,¯yi−θ⊤¯ri}. 6Chapter 1. Introduction This linear optimization problem penalizes the (training) errors ei=θ⊤¯ri−¯yi,i= 1, . . . , N in an asymmetric way. The positive errors are weighted with coefficient (1 −τ) and the negatives with coefficient τ. If τ∈(0,1) is close to zero, then the positive errors will be highly penalized (in comparison with the negative ones). This means that every optimal solution θτto the linear optimization problem will tend to make most of the errors negative. This implies that θ⊤ τrkcould be used as a probabilistic lower bound for yk. In a similar way, a probabilistic upper bound could be obtained taking τ∈(0,1) close to 1. Under rather mild assumptions, any minimizer θτof the proposed optimization problem can be used to obtain an estimation of the τquantile. That is, θ⊤ τrkserves as an estimation of the τ quantile associated with yk. See [65], [71] and [66] for further details. One of the main limitations of quantile regression is that a large number of training samples Nis required if one desires to obtain probabilistic guarantees of the method when τis chosen close to the extremes of the interval (0,1). This is due to the fact that estimating the probability of rare events requires a large number of samples. For example, the number of independent identically distributed samples required to obtain the 1 −ϵquantile of a monovariable random variable grows with 1 ϵ(see [72], [61] and [62]). Comparatively, computing multivariate prediction regions becomes a harder task. The simplest way would be to obtain intervals considering each variable independently and then construct box-shaped regions (see section 2.2.3 in [73]). 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 [74, 75] 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 a p−step prediction instead of a multivariate system. In a similar manner, [76, 75] calculate prediction regions for Vector Auto Regression (VAR) models in the field of econometrics. Following similar methodologies, [77] provides intervals for an entire path of forecasts in Markov processes. The problem of refining previously estimated prediction regions in order to attain a coverage probability closer to the desired one is researched in [78] and [79]. However, all these aforementioned techniques do not deal with multivariate output systems. Instead, the computed regions corresponds to the prediction of the same variable for different time steps. In the field of multivariate output systems, the literature is rather scarce. For example, in [80], prediction regions for a simple multivariate linear regression model are obtained. On the other hand, [81] 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 ellipsoidshaped region can be obtained. Also, [82] manages to obtain ellipsoidal regions for dynamical systems by means of an Inverse Regression (IR) scheme. This IR scheme is more efficient and reliable than the classic regression approaches when 1.1. Motivation and objectives 7 high dimensional data is available. In [83], 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 [84] are used to obtain prediction regions for any method producing a central prediction of the outputs. Similarly, it is proposed in [85] a framework to obtain a stochastic model based on a deterministic model of a dynamic system in the field of robotics. On the other hand, assuming that any finite set of samples follows a multivariate normal distribution, Gaussian processes (GPs) [23] 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. 1.1.3 Model predictive control Once reached this point, it is important to recall that the purpose of finding a good model of a dynamic system is because we want to induce a certain behaviour to these systems, that is, we want to control them. In this dissertation, Model Predictive Control (MPC) [86, 87] has been chosen as the control algorithm to be used in the proposed control examples. MPC corresponds to a set of computer control techniques sharing some common ideas. Its development exploded at the end of the 1970s thanks to the works of [88, 89], receiving attention from both academia and industry. In these works, they used predictions obtained by means of finite impulse response models or truncated step response models to obtain optimal control actions while minimizing the tracking error of the output of the system. However, most of these algorithms had an heuristic nature, they were not able to deal with disturbances and lacked stability guarantees. This first generation of predictive controllers was followed by the Generalized Predictive Control (GPC) algorithm [90, 91] which comprised these previous controllers based on input-output models. On the other hand, the state-space interpretation of MPC [92, 93] appeared, and quickly became the standard formulation in MPC. A good review of some properties of the MPC controllers, including stability can be found in [94]. The main ideas around MPC are: •MPC is an optimal controller in the sense that the input is calculated such that it minimizes a certain cost function, mainly penalizing the tracking error with respect to the desired reference and the control effort. This error can be computed thanks to a mathematical model of the system, which leads us to the following point. •As the name suggests, an accurate mathematical model describing the behaviour of the system is needed. It is possible to use whether input-output 8Chapter 1. Introduction models or state-space models. Because of the numerical nature of the controller (i.e. an optimization problem is solved to compute the value of the input instead of obtaining an explicit control law), discrete-time models are usually used. Furthermore, as it is easier to guarantee the stability of the controller using state-space models, these are also widely used. •Another interesting point of the MPC controllers is their ability to tackle constrained control problems easily. As the input is obtained by solving an optimization problem, it is very natural to add constraints. These constraints may reflect limited capability in the actuators, some physical limits of the system or even tackle security considerations. •The last idea would be the receding horizon scheme. This means that even though the value of the input is computed for many future time instants, only the first one (that is, the one corresponding to the actual time instant) is applied, discarding the rest. Then, at the next sampling step, a new whole set of inputs is computed taking into account the new available information. Note that, if we applied the whole set of values computed at the first time instant, this would be open-loop control. Thus, the receding horizon scheme provides feedback to the controller. From these properties of the MPC controllers, it is easy to see the many advantages that they present, i.e. •The control problem is formulated in the time domain, in a flexible and intuitive manner in contrast to other control formulations that need to be designed in the frequency domain. Also, it is applicable to any system without regard to its open-loop stability, delays, etc. •In the most general case, it allows to consider linear and nonlinear, univariate and multivariate systems alike using the same formulation of the controller. •It presents delay compensation explicitly and measurable disturbances can also be easily compensated. •It is possible to take into account the knowledge of the evolution of the reference in order to track the signal to the reference before it changes. •It is possible to deal easily with constraints. However, there are also some drawbacks •A sufficiently precise model of the system to be controlled is needed. That is, the performance of the controller depends heavily on the quality of the model. An inappropriate model of the system can lead to undesired performance of the controller in closed-loop operation. •It is necessary to solve an optimization problem at every time instant k. Even though nowadays the computational power of the computers have 1.1. Motivation and objectives 9 been increasing enormously, many nonlinear MPC (NMPC) or robust MPC problems still remain intractable in practice due to the high computational cost of the algorithms. As it was discussed previously, it is very common that the dynamics of the process to control are almost completely unknown, forcing us to rely on black-box models. Anyway, in most of these cases, there exists an appropriate model structure in the literature that is able to reflect the dynamics of the system and thus finding a good model is not an impossible task. On the other hand, with respect to the computational cost, there are several ways to address these problems, like developing specific optimization algorithms to solve MPC problems [95, 96] or resorting to sub-optimal MPC schemes [97, 98]. Another way would be to rely on explicit MPC implementations [99, 100]. Here, an analytical solution of the MPC controller can be obtained if the optimization problem can be casted as a certain class of quadratic programming (QP) problems [101, 102]. The obtained control law is piece-wise linear, that is, different linear state feedback control policies are obtained for different polyhedral regions. However, the number of regions grows up largely with respect to the number of states, constraints, etc. becoming intractable for many MPC problems. Besides the aforementioned classic MPC schemes, a huge number of MPC variants have been researched over time, i.e. approaches which considers disturbances and model mismatches. Actually, the classic deterministic MPC approach presents some inherent robustness [103, 104], which can be studied by means of the Input- to-State stability (ISS) [105] or robustly asymptotically stability [106] schemes. However, it is not possible to guarantee that the constraints are fulfilled in a robust manner. For that purpose, a description of the uncertainties is necessary. This led to min-max MPC, where the worst-case values of the disturbances were taken into account to compute the control actions [107, 108]. One step further was to optimize feedback control laws instead of control actions to improve the robustness of the controller [109]. The main drawback of this subclass of MPC controllers is that their online computational time is incredibly high, sometimes becoming intractable, especially when optimizing feedback control laws. Since then, numerous robust MPC schemes have been proposed. The most important one is probably the tube-based MPC [110, 111, 112]. Here, a nominal control action is computed alongside a proportional policy. As the policy is assumed to be a proportional gain, the optimization problem is not infinite dimensional and thus it is generally tractable. Also, one could rely on constraint tightening schemes [27] to fulfill the constraints under the presence of any kind of disturbance. However, robust MPC is, in general, very conservative because it addresses always the worst case, that is, it is always considering the extreme values of the disturbances. For that reason, stochastic MPC [113, 114] and scenario approaches [115, 116] which considers the probability distribution of the disturbances have 10 Chapter 1. Introduction attracted the attention of many researchers. These techniques allow us to develop controllers that fulfill the constraints with a specified probability, i.e. relaxing the worst case optimization of robust MPC. For example, in the scenario approach, this can be done by taking a finite number of realizations of the disturbance during the online computation of the control action. In this approach, however, the number of scenarios to be generated increases with the dimension of the problem, leading to unaffordable computational times in many cases. As most of these approaches complicate the online optimization problem to be solved online, one could wonder if it is possible to obtain a stochastic MPC controller while keeping the online optimization problem as simple as the deterministic MPC problem. This leads us to probabilistic validation approaches [117, 118] where, by means of offline simulations of the closed-loop system, it is possible to determine if the controller fulfills some probabilistic specifications. This means that the computational burden is outside the control loop and thus the online optimization problem can be solved easily. Also, it can be applied to any controller, even if it is not an MPC controller. However, this method only provides a “yes” or “no” answer without assessing numerically the performance of the controller, which makes it hard to compare with other controllers. In some cases, it could even lead to non feasible solutions if none of the controllers satisfies the probabilistic constraint. Other popular variants of the standard MPC controller are the following: •Economic MPC [119, 120] allows to optimize different cost functions in order to improve the profitability of a certain process. The stability results of standard MPC are not, in general, applicable for economic MPC. However, under some assumptions [121], it was proven that it is possible to find a Lyapunov function of the closed-loop system in order to ensure stability. •Distributed MPC [122, 123] considers a decentralized control system, reducing the computational burden of the individual optimization problems and becoming more scalable and robust to failures, which may be helpful when controlling large-scale systems. However, the performance may be worse than the performance of the centralized controller. •MPC for tracking [124, 125] tackles the problem of arbitrarily changing the reference signal of the controller. This approach preserves stability and guarantees that the system can be steered to any admissible equilibrium point, no matter what the initial equilibrium point is, without losing the recursive feasibility. This strategy also provides a larger domain of attraction for a given prediction horizon and, in the case of unreachable references, asymptotic convergence to the closest reachable reference is proven as well, only at the expense of adding a few more decision variables. 1.2. Outline 11 1.1.4 Objectives of this dissertation Many different research fields have been presented throughout this introductory chapter: system identification, interval predictions, MPC, etc. This dissertation is oriented towards developing new methods that may improve the performance of the state-of-the-art methods in each one of the previously presented fields. Summarizing, the objectives of this dissertation are: 1. Developing new forecasting schemes for time-series and nonlinear systems. It is interesting to obtain not only the expected output but also a region where the output can be found with a certain probability. 2. Proposing new data-driven black-box models to be used in MPC strategies. As the proposed models may reflect better the behaviour of certain systems, they might increase the performance of certain controllers. 3. Obtaining probabilistic bounds on the constraint violation level for any kind of controller in a stochastic setting. These bounds could be used to tune the hyper-parameters of a controller in a soft-constraint approach. Taking into account these objectives, the following results were obtained: 1. Corresponding to the first objective, some predictors based on dissimilarity functions are proposed. These include nominal predictions, interval predictions and prediction regions. All these methods were tested against some baseline techniques in many examples and showed a better performance. 2. For the second objective, it is shown that, by means of the aforementioned dissimilarity function, it is possible to identify a model from the available input-output data. Through different numerical examples, it can be seen that the model performs better than other machine learning techniques. The model has been used within a MPC for tracking framework to show the effectiveness of the approach in control problems. 3. Finally, two different sharp bounds on the constraint violation rate were obtained. These bounds were used in the context of data center energy optimization to quantify the quality of service provided to the users, proving its usefulness. 1.2 Outline All the results of this dissertation converge to MPC at some point, that is, they are useful to improve the performance of MPC controllers or they allow to develop new MPC schemes that may attain better results than other state-of-the-art MPC controllers. For the sake of simplicity, the contributions are divided into three groups: novel algorithms for prediction and forecasting, new modeling methods 12 Chapter 1. Introduction for controller design and, finally, randomized algorithms for stochastic control. Having this idea in mind, the text was divided into three different parts. Part I presents the developed predictors based on dissimilarity functions. These proposed methodologies are related to the ones appearing in the context of Direct Weight Optimization [57] and the Kriging method [126, 127]. First, the concept of dissimilarity function is exposed, along with the proposed dissimilarity function for this dissertation. From that, a nominal predictor can be easily obtained. Given a regressor, a prediction is obtained as a combination of past outputs of the system using some weights that are obtained by means of an optimization problem. After that, a method to tackle the problem of uncertainty quantification is proposed. Again, by means of the aforementioned dissimilarity functions, it is shown that it is possible to obtain interval predictors for the univariate case and prediction regions for the multivariate case. Part I is divided into two chapters: Chapter 2 works as an introduction to chapter 3, introducing the basics of dissimilarity functions and proposing a first predictor to obtain the expected value of a time-series. A numerical example of the forecasting of the Dow Jones Industrial Average Index is shown to verify the good performance of the algorithm. Chapter 3 presents the methodologies to quantify the uncertainty of the predictions. First, a method to obtain an empirical probability density function for a univariate function is discussed. From this empirical probability distribution, it is possible to compute quantiles of a desired probability, obtaining an interval predictor. Later, the methodology is extended to tackle the multivariate case. As obtaining the probability distribution of a multivariate random variable is generally intractable, an alternative method bypassing the probability density function allows us to obtain prediction regions of such variables. Some numerical examples are presented in order to show the improvements with respect to other baseline techniques. On the other hand, part II presents novel methodologies to obtain models of nonlinear systems or time-series. Here, the methodologies proposed in part I are extended in order to obtain a model of the system by feed backing the 1-step ahead predictions and regularizing the optimization problem used to compute such predictions. Two different models are proposed in this chapter: a linear time varying (LTV) model obtained by weighting appropriately the local data within the data set and a kernel-based version of the state-space kriging where the non-linearity is modeled by means of a kernel function. The obtained models can be used easily to make predictions or design MPC controllers. Forecasting and Control examples are both provided to show the effectiveness of the proposed approaches. Part II is divided into two chapters: 1.3. Publications 13 Chapter 4 presents the state-space kriging approach, named after the kriging method due to its similarities and also because the weights used to obtain the predictions become the state of the new model. This chapter is focused on autonomous systems and time-series. Thus, it presents methodologies to obtain a model of the outputs from past data of the system that can be used to make predictions. Numerical examples with comparisons are provided in order to show the effectiveness of the proposed approaches. Chapter 5 introduces the state-space kriging for non-autonomous systems. It is shown that the input can be considered easily within the optimization problem with very few changes. As the approach considered in this chapter can tackle the problem of forecasting systems with manipulable inputs, the proposed approach is suitable to design, for example, MPC controllers whose prediction model is a state-space kriging model. The stability of the proposed controller is proven and a numerical example and a real experiment are provided to verify the performance of the controller. Finally, part III presents novel bounds on the constraint violation rate that can be used to quantify the performance of different controllers in a stochastic setting in order to choose the most appropriate controller among them. These controllers are not necessarily MPC controllers, that is, any controller can be assessed by means of these techniques. Another advantage of these methods is that the computational burden of the process is completely offline and thus the online computation time of the controller does not change. The proposed methods are finally validated in the context of the management of a data center. Part III is divided into two chapters: Chapter 6 presents the concepts needed to obtain the bounds on the empirical constraint violation level based on the results by Chernoff in [128] and by Alamo et al. in [61]. These results provides two different expressions that can be used to measure the reliability of any controller. Chapter 7 introduces the model of a data center (dealing with both the thermal modeling and the task model for the computers) that will be used as a benchmark for the previously developed bounds. In this chapter, not only the model of the data center and the implied constraints are presented, but a management approach based on MPC is also proposed. As the problem for this specific system is generally hard to solve (integer optimization, non-convex model, etc.) the solution of the problem is computed by means of a particle-based solver implemented in a Graphic Processing Unit (GPU). 1.3 Publications The results shown throughout this thesis have been published in several journals and conference proceedings, although some of them are still under review. 20 Chapter 2. Forecasting using dissimilarity functions coordinate system. However, many of the dissimilarity functions that can be found in the literature are not invariant. For example, any dissimilarity function based on the distance of dto the elements of D, such as that of equation (2.1), will be dependent on the particular choice of the coordinate system. The proposed optimization problem (2.2) is a strict convex optimization problem subject to convex constraints. This means that it has a unique solution [139]. Note that the numerical resolution can be addressed using a dual formulation. In the dual formulation for this particular optimization problem, the number of dual decision variables is equal to the number of equality constraints (n+ 1) which is in many situations much smaller than the number of primal variables (N). On the other hand, the gradient of the objective function in the dual formulation can be obtained in a direct way because once the dual variables are fixed, the optimal values for the primal variables are obtained by solving Nseparable optimization problems. The computations in this dissertation have been made using an accelerated gradient method in the dual variables (see [140], [141] and [142]). As it is formally stated in the following property, the optimization problem has an explicit solution for the particular case γ= 0. Property 2.7. J0(d, D) has the following explicit expression J0(d, D) = N−1+ (d−dc)⊤(DD⊤−Ndcd⊤ c)−1(d−dc), where dc=N−1D1⊤and 1⊤is a column vector with all its Ncomponents equal to 1. Proof. The optimization problem is solved using a dual formulation where µ∈Rn denotes the multipliers associated with the equality constraint d= N X i=1 λ(i)di=Dλ, and νis the multiplier corresponding to the equality 1 = N X i=1 λ(i)=1λ. the Lagrange function is L(λ, µ, ν) = λ⊤λ+µ⊤(Dλ −d) + ν(1λ−1). Denote λ∗,µ∗and ν∗the optimal values for the primal and dual variables. From ∂L(λ∗,µ∗,ν∗) ∂λ = 0, it is clear that the optimal vector λ∗is given by λ∗=−1 2(D⊤µ∗+1⊤ν∗).(2.3) 2.1. Proposed dissimilarity function 21 Since 1λ∗= 1, D1⊤=Ndcand 11⊤=N, it is possible to premultiply both terms of the last equality by 1to obtain 1 = −1 2(1D⊤µ∗+Nν∗) =−N 2(d⊤ cµ∗+ν∗). Therefore, ν∗=−2 N−d⊤ cµ∗. Substituting the expression for ν∗in (2.3) yields λ∗=−1 2D⊤µ∗−1⊤(2 N+d⊤ cµ∗) =1⊤ N−1 2(D⊤−1⊤d⊤ c)µ∗.(2.4) Premultiplying by D, the following equation is obtained Dλ∗=dc−1 2(DD⊤−Ndcd⊤ c)µ∗.(2.5) From the equality constraint Dλ∗=dand (2.5), the value of µ∗can be obtained as µ∗=−2(DD⊤−Ndcd⊤ c)−1(d−dc). Substituting µ∗in equation (2.4), it is possible to compute λ∗as λ∗=1⊤ N+ (D⊤−1⊤d⊤ c)(DD⊤−Ndcd⊤ c)−1(d−dc). Finally, taking into account that 1(D⊤−1⊤d⊤ c) = (Nd⊤ c−Nd⊤ c)=0, the following expression is obtained (D⊤−1⊤d⊤ c)⊤(D⊤−1⊤d⊤ c) = DD⊤−Ndcd⊤ c. From last equality and the expression for λ∗, the value of J0(d, D) can be computed as J0(d, D)=(λ∗)⊤λ∗=N−1+ (d−dc)⊤(DD⊤−Ndcd⊤ c)−1(d−dc).■ The previous result shows that the dissimilarity function is a quadratic function on the argument dfor the particular case γ= 0. For the more general case in which γ > 0, it is possible to infer from the Karush-Kuhn-Tucker (KKT) optimality conditions [139] that the dissimilarity function Jγ(d, D) is a piecewise convex quadratic function with respect to d. 22 Chapter 2. Forecasting using dissimilarity functions 2.1.1 Clarifying example Figure 2.1 shows a cloud of points obtained from a certain probability distribution. On the left, the original cloud of points is shown and, on the right, this cloud of points is rotated an angle π/4. Red markers correspond to the points in the data set Dwhereas blue markers correspond to two given values of d, i.e. d1=0.25 1, d2=0.75 1. Note that these values correspond to the left figure, where they are not rotated. Assuming that γ= 0, the obtained values of the dissimilarity for these vectors d1 and d2are J0(d1, D)=0.0056, J0(d2, D)=0.0191. It is easy to see that the dissimilarity is smaller for d1because the concentration of data points is larger in that region, in contrast to d2where the concentration of red markers is small. This means that d1is a more likely event than d2. -0.5 0 0.5 1 1.5 -0.5 0 0.5 1 1.5 2 2.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Figure 2.1: Left: original cloud of points. Right: rotated cloud of points. Also, note that, as it was proven in property 2.6, the value of the dissimilarity function is invariant with respect to affine transformations. 2.2 Dissimilarity functions and regression This section shows how dissimilarity functions can be used in the context of regression. Denoting the regressor as z∈Rnzand the output as y∈Rny, it is assumed that the following data sets Zand Yare available Z=¯z1¯z2. . . ¯zN, Y =¯y1¯y2. . . ¯yN, where ¯yi∈Rnyand ¯zi∈Rnz∀i= 1, ..., N are past samples of the outputs and the regressor respectively. 2.2. Dissimilarity functions and regression 23 Fixing a certain zand computing Jγ z y,Z Ygives us the dissimilarity of a certain output yconditioned to the value of the regressor z. Thus, the dissimilarity function becomes Jγ z y,Z Y= min λ(1),...,λ(N) (1 −γ) N X i=1 λ2 (i)+γ N X i=1 |λ(i)|(2.6a) s.t. z= N X i=1 λ(i)¯zi(2.6b) y= N X i=1 λ(i)¯yi(2.6c) 1 = N X i=1 λ(i).(2.6d) Given z, consider, from all the possible values of y, the one that minimizes the dissimilarity function defined by the previous optimization problem ˜y(z) = arg min yJγ z y,Z Y.(2.7) Note that minimizing the dissimilarity corresponds to maximizing the similarity. Thus, ˜y(z) is the most likely value to happen given the regressor zand the data sets Zand Y. This ˜y(z) could be used as a forecasting of the next output. Now, instead of Jγ z y,Z Y, consider Jγ(z, Z), that is, Jγ(z, Z) = min λ(1),...,λ(N) (1 −γ) N X i=1 λ2 (i)+γ N X i=1 |λ(i)|(2.8a) s.t. z= N X i=1 λ(i)¯zi(2.8b) 1 = N X i=1 λ(i).(2.8c) As the optimization problem in (2.8) is a less constrained version of that in (2.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. 24 Chapter 2. Forecasting using dissimilarity functions Now, denote the vector of weights λthat minimizes the optimization problem in (2.8) as λ∗. From this optimum vector of weights λ∗, ˜y(z) can be computed as ˜y(z) = Y λ∗= N X i=1 λ∗ (i)¯yi.(2.9) Then, it is easy to see that Jγ(z, Z) = Jγ z ˜y(z),Z Y and thus ˜y(z), computed using the optimal weights from (2.8), is the one that minimizes the dissimilarity function for a fixed z(see equation (2.7)). Note that this forecasting is made by means of a linear combination of past outputs whose weights are computed obtaining an affine envelope of the regressor z. As it is explained 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). It is easy to see that when every λ(i)≥0,∀i= 1, ..., N, then PN i=1 |λ(i)|= 1. However, when some of the components of λare negative, PN i=1 |λ(i)|>1. In other words, the cost term PN i=1 |λ(i)|becomes larger when extrapolating points (i.e. using negative values of λ(i)). This means that convex combinations of λ(i)are encouraged and thus interpolation is preferred, which may improve the predictions when nonlinear systems are taken into account. The optimal value for γ≥0 could be obtained in different ways. A reasonable choice could be, for example, to choose the one that minimizes the error of the predictions with respect to the real outputs in a validation set. As it is stated in the following property, the estimation obtained for the particular case γ= 0 matches the one given by the linear least squares method. Property 2.8. Given z, the estimation ˜y(z) = arg min yJ0 z y,Z Y, matches the estimation obtained by linear least squares. Proof. From equation (2.9), the optimal value for the estimation is ˜y(z) = ΣN i=1λ∗ (i)¯yi, where for the particular case γ= 0, λ∗ (i),i= 1, . . . , N, are the optimal values of 2.2. Dissimilarity functions and regression 25 the optimization problem min λ(1),...,λ(N) N X i=1 λ2 (i) s.t. z= N X i=1 λ(i)¯zi 1 = N X i=1 λ(i). Defining R=¯z1¯z2. . . ¯zN 1 1 . . . 1, r =z 1, the equality constraints can be rewritten as Rλ =r. (2.10) From the KKT optimality conditions, the optimal solution is given by (see subsection 10.1.1 in [139]) INR⊤ R0(nz+1)  λ∗ φ∗=0 r, where φ∗corresponds to the optimal dual decision variables corresponding to the equality constraint (2.10). The previous equation can be rewritten as λ∗=−R⊤φ∗ Rλ∗=r. From here, it is inferred that −RR⊤φ∗=r, which implies φ∗=−(RR⊤)−1rand thus finally obtaining λ∗=R⊤(RR⊤)−1r. Therefore, ˜y(z) = Y⊤λ∗=Y⊤R⊤(RR⊤)−1r=r⊤(RR⊤)−1RY. This corresponds to the least squares estimation obtained when the considered regressors are the vectors ¯z⊤ i1⊤,i= 1, . . . , N (see [143], [64]). ■ Previous property shows that the proposed estimation method encompasses the least squares method for the particular case γ= 0. A family of optimal estimators is obtained if γis considered a tuning parameter. From now on, in order to make the manuscript more readable, the notation will be simplified by removing Zand Yfrom the dissimilarity function (since they are assumed to be fixed). That is, the notation becomes Jγ z y,Z Y=Jγ(z, y). 26 Chapter 2. Forecasting using dissimilarity functions Moreover, the minimal value of the dissimilarity function given zwill be denoted as J∗ γ(z), that is J∗ γ(z) = Jγ(z, ˜y(z)) = min yJγ(z, y). The following remark summarizes the process to use the dissimilarity functions in the context of regression. Remark 2.9. Given the regressor z, the dissimilarity function Jγ(·,·)and the data sets Zand Y, solve optimization problem in (2.8) to obtain λ∗. Once λ∗has been obtained, the forecasting ˜y(z)is computed by means of equation (2.9). 2.3 Application: forecasting stock prices using dissimilarity functions In this section, a numerical example is presented to show the effectiveness of the proposed dissimilarity function as a predictor. Consider the evolution of the price of a stock as a time series pk∈ P, where k is a time unit and Pis the range of values that the stock price can take. The state of this time series can be described as the value at time kof a series of technical indicators. These technical indicators can be past values of the price, stock price returns or more complex metrics in general. Thus, this group of technical indicators will form the regressor zk. The objective is to be able to predict up to l-steps (i.e. l-days) ahead the price of a stock, that is, to obtain ˜pk+lat time k. Thus, ykcorresponds to yk=pk+1 pk+2 . . . pk+l⊤. As it was shown in the previous section, Zand Yare compounded of past samples of zand y, denoted as ¯ziand ¯yi. The approach considered in this example uses only a small subset of the database, denoted as Ω(zk), to compute the predictions ˜pk+l. This subset Ω(zk) is chosen by finding the NΩclosest points to zkwithin the data set by means of a certain distance measurement. Thus, Ω(zk) is compounded of NΩpairs (¯zi,¯yi) of the full data set. Once the data to be included in Ω(zk) are selected, the optimization problem (2.8) is solved with this reduced data set. Also, note that the forecasting of pk+lis also made with this small portion of the complete data set. Algorithm 4 gives a formal description of the proposed approach. Remark 2.10. The distance can be any measure of how close are zkto ¯zi. Euclidean distance would be the most typical choice, but it is also possible to consider other aspects like the time span between samples. In this way, recent data could 2.3. Application: forecasting stock prices using dissimilarity functions 27 Algorithm 1: l-step ahead stock forecasting using local data Data: Z,Y,zk,NΩand γ. Result: ˜pk+1 ˜pk+2 . . . ˜pk+l. 1Compute the distance for each ¯ziin the database to zk; 2Create a list of the entries in the data base sorted according to the computed distances. Denote ¯zjand ¯yjas the regressor ¯ziand the output ¯yiof the j-th entry in this ordered list respectively; 3Build Ω(zk) using the first NΩentries of the ordered list, that is, Ω(zk)≜{(¯zj,¯yj)} ∀j= 1, . . . , NΩ. 4Solve the following optimization problem: min λ(1), ..., λ(NΩ) (1 −γ) NΩ X j=1 λ2 (j)+γ NΩ X j=1 |λ(j)| s.t. NΩ X j=1 λ(j)= 1, NΩ X j=1 λ(j)¯zj=zk. and compute the forecasted prices as: ˜pk+1 ˜pk+2 . . . ˜pk+l⊤= ˜yk= NΩ X j=1 λ(j)¯yj. be prioritized when selecting the elements of Ω(zk). Other aspects like seasonality could also be taken into account. 2.3.1 Results The data set to be used in this example was obtained from the data provider Bloomberg and it is composed of the daily closing prices of the Dow Jones Index from 2005 to 2016. The data was divided into a training set (Zand Y) from 2005 to 2014 and a test set from 2015 to mid-2016. This period was chosen because the market does not follow a certain trend (bullish or bearish) that would make the forecasting trivial. In order to reduce the noise of the time series, the prices are smoothed using a 5-day Exponential Moving Average (EMA) computed as: pm EMA,k =2 m+ 1pk+1−2 m+ 1pm EMA,k−1, 28 Chapter 2. Forecasting using dissimilarity functions with pm EMA,0=p0and m= 5 in this case. Note that the smoothing applied here is very light in comparison to the usual values of mused in short-term forecasting (12 and 26 day EMA [144]). This preserves fast price fluctuations at the cost of making the forecasting process harder. The regressor zkis composed of the last ten days smoothed prices, as well as the 5-day and 10-day relative difference percentage of unsmoothed prices (RDP) [145] i.e. RDPm k= 100 pk−pk−m pk, being mequal to 5 and 10 respectively. Thus, zk=pkpk−1. . . pk−9RDP5 kRDP10 k⊤. The size of Ω(zk) is set to NΩ= 250 and γ= 0. The forecast and real prices can be seen in figures 2.2 and 2.3. The results show that the forecasting is accurate enough at first and it becomes worse as the prediction horizon increases, as it can be expected. In order to compare the results obtained with the proposed approach, the results obtained with a persistence predictor and a Neural Network (NN) are shown. This persistence predictor works as a martingale, that is, ˜pk+l=pkwhereas the NN corresponds to a multi-layer perceptron with 20 neurons in the hidden layer and trained following the Levenberg-Marquardt rule. The numerical results are shown in table 2.1. In particular, the mean squared errors (MSE) of the proposed approach and the aforementioned baselines are shown in table 2.1a. It can be seen that the proposed approach achieves the lowest errors. However, note that for long step ahead predictions, the persistence predictor achieves results as good as the other techniques, a consequence from the long term random walk nature of the financial market. Table 2.1: Results of the Dow Jones forecasting example. (a) Mean squared errors (MSE). kProposed MLP Persistence 1 3,294.8 3,577.0 5,296.4 2 12,433.2 14,190.8 16,867.4 3 26,659.0 30,829.6 31,975.0 4 45,475.9 50,794.9 49,072.4 5 66,859.6 77,017.5 66,943.4 (b) Standard deviation of the errors (σ). Proposed MLP Persistence 57.4 59.9 72.8 111.5 119.3 130.0 163.1 175.6 178.9 212.9 224.9 221.6 257.8 276.0 258.8 On the other hand, table 2.1b shows the standard deviation of both the proposed approach and the baselines. Again, the proposed approach attains the lowest standard deviation in the errors, which means that the errors are not only smaller in mean but also they are more concentrated. 2.3. Application: forecasting stock prices using dissimilarity functions 29 0 50 100 150 200 250 300 Day 1.55 1.6 1.65 1.7 1.75 1.8 1.85 Closing Price 104 Forecast Real 0 50 100 150 200 250 300 Day 1.55 1.6 1.65 1.7 1.75 1.8 1.85 Closing Price 104 Forecast Real 0 50 100 150 200 250 300 Day 1.55 1.6 1.65 1.7 1.75 1.8 1.85 Closing Price 104 Forecast Real Figure 2.2: Forecasted and real prices (5-day EMA) for 1 to 3 days forecasting. 36 Chapter 3. Probabilistic prediction regions Algorithm 2: Interval estimation [y− τ(z, γ, c), y+ τ(z, γ, c)]. Data: z,τ∈(0,1), γ≥0, c≥0, Z,Y,ˇ Y. Result: y− τ,y+ τ. 1Obtain the dissimilarity function (see Definition 2.1) for each element of ˇ Y: aj=Jγ(z, ˇyj), j = 1, . . . , M. 2Compute the conditioned probabilities (see equation 3.1): pj= ecpγ,c (z, ˇyj) = e−caj M P ℓ=1 e−caℓ , j = 1, . . . , M. 3Compute the indexes ℓ+ τand ℓ− τcorresponding to the lower and upper conditioned τ-quantiles (see (3.2) and (3.3)): ℓ+ τ= smallest integer ℓsatisfying ℓ P j=1 pj≥1−τ, ℓ− τ= largest integer ℓsatisfying M P j=ℓ pj≥1−τ. 4Make y− τ= ˇyℓ− τand y+ τ= ˇyℓ+ τ. Consider now the role of parameter c≥0 in the discrete empirical conditioned distribution given in equation (3.1). On the one hand, the choice c= 0 provides a flat distribution in which each element of ˇ Yhas a conditioned probability equal to 1 M. On the other hand, large values of cprovide narrow distributions centered around the point in ˇ Ythat minimizes, given z, the dissimilarity function Jγ(·,·), i.e. ˜y(z). Consequently, for a fixed value of γ, larger values of creduce the size of the obtained interval at the expense of increasing the fraction of outputs that are not contained in the interval estimations. Therefore, given γ, the corresponding value for cshould be chosen as the largest value of cthat guarantees in the validation set that the obtained intervals contain the outputs with the desired probability. From the discussion above, it is clear that the parameter ccorresponding to a particular choice of γ > 0 (denoted cγ) is determined by τ. As it is detailed in Algorithm 3, cγis chosen as the largest value of c(up to a given accuracy ϵ > 0) that guarantees in the validation set that the obtained confidence intervals contain the outputs with the desired probability, that is, no smaller than 1 −2τ. Parameter γ > 0 can be obtained by maximizing the likelihood ratio which, for a 3.1. Univariate case 37 Algorithm 3: Optimal value of c≥0, for given γ≥0 and τ∈(0,1) Data: τ≥0, γ≥0, cmax >0 and ϵ > 0, Z,Y,ˇ Yand the validation set V. Result: cγ. 1cmin = 0; 2while cmax −cmin ≥ϵdo 3c=1 2(cmax +cmin); 4Compute, using Algorithm 2, the NVintervals Ii= [y− τ(¯zi, γ, c), y+ τ(¯zi, γ, c)], i = 1, . . . , NV. 5Make n+ viol equal to the number of violations of the upper constraints ¯yi≤y+ τ(¯zi, γ, c), i = 1, . . . , NV, and n− viol equal to the number of violations of the lower constraints ¯yi≥y− τ(¯zi, γ, c), i = 1, . . . , NV. 6if max{n+ viol, n− viol} NV < τ then 7cmin =c; 8else 9cmax =c; 10 end if 11 end while 12 cγ=cmin; 38 Chapter 3. Probabilistic prediction regions specific γand corresponding cγ, is defined as Lγ= NV X i=1 log ecpγ,c (¯zi,¯yi). Thus, using ˇ Y={ˇy1,...,ˇyM}, the optimal value of γis given by γ∗ τ≈arg max γ∈Γ NV X i=1 log      e−cγJγ(¯zi,¯yi) M P j=1 e−cγJγ(¯zi,ˇyj)      ,(3.4) where Γ is a set containing all the possible values considered for γ. Remark 3.2. Other criteria can be used to compute γ∗ τ. For example, γ∗ τcould be obtained by minimizing a cost function penalizing the average length of the intervals and/or the average prediction error with respect to the conditioned median introduced in Remark 3.1. However, explicitly minimizing the size of the intervals may translate into an increased violation rate when the validation set has not a sufficiently large number of samples, hence not being representative enough of the real distribution of y. 3.1.2 Clarifying example: uniform distribution A sample of 600 points in Ris obtained from a uniform probability function with support [0,1]. One half of the available points is used as a training set and the other half is used as a test set. Figure 3.1 shows the empirical probability density functions estimated using different values of the parameters cand γ. In this case, c= 1.5 and γ= 5 is the pair that achieves the best fit for the distribution proposed in this example. 3.1.3 Numerical example: Lorenz attractor The Lorenz attractor is a system of ODEs known for having chaotic solutions with certain values of the parameters of the system. The equations that define the system are the following do dt=σ(p−o) dp dt=o(ρ−q)−p(3.5) dq dt=op −βq , where σ,ρand βare real scalar parameters. In this example, these parameters take the values σ= 10, ρ= 28 and β= 8/3. Furthermore, in order to obtain the 3.1. Univariate case 39 -0.5 0 0.5 1 1.5 c=150, =0 0 0.5 1 1.5 -0.5 0 0.5 1 1.5 c=15, =0.4 0 0.5 1 1.5 -0.5 0 0.5 1 1.5 c=7.5, =1 0 0.5 1 1.5 -0.5 0 0.5 1 1.5 c=1.5, =5 0 0.5 1 1.5 Figure 3.1: Estimated probability distribution functions. necessary data, the ODEs have been integrated numerically with a fixed time step of Ts= 0.1sand initial conditions o0= 1, p0= 1 and q0= 1. Here, it is considered the task of forecasting the one-step ahead value of o, i.e., yk+1 =ok+1, using the two previous values of o, that is, the regressor vector will be zk= [yk, yk−1]⊤. To start with, 2500 data points are considered, normalized in the [0,1] range. Different sizes for the data sets Zand Yare considered in this example (200, 350 and 500 points). The validation set Vconsists of 1000 data points and other 1000 data points are used as a test set, denoted by T. Here, we consider the dissimilarity function of remark 2.3. The set Γ is taken from [0,3] using a 0.1 sampling step. On the other hand, ˇ Yis obtained from a grid of equally distant points in the interval [−0.1893,1.2298] sampled with a 1.4191 ×10−4step. Two different techniques will be considered as benchmarks. The first one is quantile regression [65], [66], a classical method for the estimation of conditioned quantiles. The second one is the set-membership method described in [49, 50]. This technique is a well-known method to generate interval bounds for a time series (usually produced from a dynamical system). For the sake of comparison, to guarantee that these bounds contain the output within a prescribed probability, the parameters ϵ, γ of [49] are chosen so that the resulting empirical probability of containing a sample within the validation set Vis no smaller than 1 −2τ. The numerical results of the proposed approach and the two benchmark techniques are shown in table 3.1 for a [o5%, o95%] interval, that is, τ= 0.05, and in table 3.2 for [o10%, o90%] (τ= 0.1). The output of the test data should be contained in the first interval with a probability of 0.9 (0.8 for the second interval). The optimal value for γhas been chosen by maximizing the likelihood function Lγ(see equation (3.4) and figure 3.3). 40 Chapter 3. Probabilistic prediction regions Table 3.1: Results for the Lorenz Attractor, interval [o5%, o95%]. Proposed QR SM NE. P. I. W. E. P. I. W. E. P. I. W. 200 0.9140 2.0578 0.8290 3.0965 0.8960 2.9378 350 0.8990 1.9352 0.8260 3.0550 0.9100 2.4773 500 0.9070 2.0223 0.8410 3.2450 0.9120 2.5671 Table 3.2: Results for the Lorenz Attractor, interval [o10%, o90%]. Proposed QR SM NE. P. I. W. E. P. I. W. E. P. I. W. 200 0.8060 1.6053 0.7450 2.2776 0.8160 2.4248 350 0.8060 1.6164 0.7270 2.0607 0.7900 1.9797 500 0.8100 1.6195 0.7630 2.4371 0.8100 2.0021 The empirical probability in the case of the quantile regression clearly does not meet the probabilistic specifications. In the case of the proposed approach and the set-membership method, the observed fraction of outputs that fall into the predicted intervals is much closer to the desired one. Note that, for all techniques, the obtained empirical probability can be below the desired probability. This could be solved relying on a probabilistic scaling scheme [63] or probabilistic validation schemes [133], [117]. Regarding the interval width, the proposed approach clearly manages to obtain the smallest intervals for each data set. For the [o5%, o95%] case, in comparison to the set-membership method, the mean interval width is 24.35% smaller. Also, it is 35.96% smaller with respect to the mean interval width obtained by means of quantile regression. On the other hand, for the [o10%, o90%] case, the intervals obtained with the proposed technique are 23.75% smaller than those obtained by means of the set-membership method and 28.21% smaller than those obtained using quantile regression techniques. Taking into account the empirical probability values and the interval widths, it is possible to conclude that the proposed approach obtains the best results. Finally, in figure 3.2, a fraction of the test set Tis shown along with the computed intervals [o5%, o95%] of the proposed approach. Note that the intervals are wider when there are trend changes in the output. Furthermore, figure 3.3 shows an example of the value of the maximum likelihood ratio Lγas a function of γ(in this case for the data set of 200 points and interval [o5%, o95%]). 3.1.4 Numerical example: Dow Jones industrial average index The example presented in the previous chapter will be used again to show the effectiveness of the proposed approach for interval forecasting. The parameters of 3.1. Univariate case 41 Figure 3.2: Test set and computed intervals for Lorenz Attractor (interval [o5%, o95%]). 0 0.5 1 1.5 2 2.5 3 -6350 -6300 -6250 -6200 -6150 -6100 L Figure 3.3: Maximum likelihood ratio as a function of γ. the algorithms are NΩ= 250, M= 1000, ˇy1= 6684.3, ˇyM= 19445, cmax = 15. As in the previous example, the dissimilarity function used is the one corresponding to remark 2.3. The set Γ is taken from [0,5] using a 0.5 sampling step. The regressor is compounded as in chapter 2. The 10-th and 90-th percentiles (i.e. τ= 0.1) corresponds to the desired probability for the price intervals. Note that, in order to compute the intervals with the proposed methodology, it is needed to consider each l-step prediction independently, obtaining ldifferent predictors, each one for a certain l-step forecasting. The results are shown in figures 3.4 and 3.5. There, the price intervals are presented as envelopes. Although these intervals are tight for l= 1, they become larger as lincreases. This is straightforward due to the fact that the uncertainty increases as the prediction horizon does. Also, note that sometimes the real price is not included within the computed price interval. This is congruent to the fact that the probability of belonging to the interval is 80%. 42 Chapter 3. Probabilistic prediction regions Figure 3.4: Price intervals for 1 to 3 days (5-day EMA). 3.1. Univariate case 43 Figure 3.5: Price intervals for 4 and 5 days (5-day EMA). 44 Chapter 3. Probabilistic prediction regions As a baseline approach, a quantile regression approach is chosen to validate the obtained results. In table 3.3 it is shown the empirical probabilities and interval width of the proposed approach and the quantile regression scheme. Although the proposed approach contains the real price with a higher probability than the desired one, the quantile regression fails to obtain intervals meeting the specified probability for any l. Table 3.3: Empirical probability and average interval width. Empirical probability Average interval width kProposed Q. regression Proposed Q. regression 1 0.8679 0.6006 162.4802 99.8091 2 0.8459 0.6101 315.7218 190.5443 3 0.8553 0.6321 477.8004 287.5774 4 0.8648 0.6761 638.8747 397.3780 5 0.8805 0.6950 813.5661 489.7816 3.2 Multivariate case In this section, the problem of multidimensional systems or time-series is tackled, providing prediction regions for the multidimensional output. As the scheme proposed in the previous section for univariate outputs is based on a numerical integration in a unidimensional space, the generalization to the multidimensional setting is not trivial because of the well-known complexity of high-dimensional numerical integration. In what follows, it is shown how to obtain prediction regions in multidimensional spaces without resorting to numerical integration, providing a fundamental advantage with respect to the previous results. First, implicit regions computed by means of the dissimilarity function are proposed. After that, an ellipsoidal approximation is proposed. Both approaches are tested in a numerical example. 3.2.1 Implicit regions In order to characterize the desired prediction region, two issues have to be tackled: the choice of the region center and the computation of the region itself. Choosing the center In the previous chapter, ˜y(z) was denoted as the value of ythat minimizes the dissimilarity function given z(see (2.7)). As a consequence of this, ˜y(z) can be considered as a good option to define the center of the region that will be obtained as it is the most likely value for the output given z,Zand Y. This optimal value can be obtained as ˜y(z) = Y λ∗where λ∗corresponds to the 3.2. Multivariate case 45 optimal solution of the strictly convex problem (2.8), that is Jγ(z, Z) = min λ(1),...,λ(N) (1 −γ) N X i=1 λ2 (i)+γ N X i=1 |λ(i)| s.t. z= N X i=1 λ(i)¯zi 1 = N X i=1 λ(i). Computing the regions In the previous section, it was shown that the dissimilarity function can be seen as a sort of surrogate of the probability distribution of ygiven z. Thus, this probability distribution peaks at ˜y(z) and decreases as the dissimilarity function increases. For that reason, the proposed prediction regions are defined as those points for which the dissimilarity function does not exceed more than a given factor αwith respect to the value corresponding to the central prediction ˜y(z), that is, J∗ γ(z). This is formally stated in the following definition. Definition 3.3 (Prediction region).For a given z,γ, data sets Z,Yand a tunable parameter α > 1, the proposed prediction region is defined as the set ∆(z) = y:Jγ(z, y)≤αJ∗ γ(z),(3.7) that is, the points ythat obtain a dissimilarity less or equal to αJ∗ γ(z). 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 problems. 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 (3.7) holds. Tuning the hyperparameters αand γ The value of γ≥0 could be obtained in such a way that the prediction error corresponding to the central prediction in a validation set is minimized. For a given γ, smaller values of αmake the regions smaller. Our objective is to obtain the smallest possible region that guarantees that a point ¯yjtaken from a validation set Vis contained in the computed region ∆(zj) with a pre-specified probability 1−τ. Again, as in the univariate case, it is assumed that the validation set V represents the true probability distribution. Then, αis chosen so that it fulfills ProbV(y∈∆(z)) ≥1−τ, 52 Chapter 3. Probabilistic prediction regions Table 3.4: Area of the regions and empirical probabilities (E.P.) 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 NArea 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 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. E.P. ≥1−τ). This means that the proposed approach provides regions of smaller uncertainty in comparison to the baselines considered. 3.3 Conclusions This chapter presented methods to obtain interval predictions and prediction regions, which can be considered as an extension of the previous chapter where only the expected value of the output of a time-series or a dynamic system was obtained. First, a method based on obtaining an empirical conditioned probability density function was proposed. However, this method turned out to be computationally expensive due to the fact that it required the numerical integration of this probability density function for every time instant. Two numerical examples were proposed: the forecasting of the closing price of the Dow Jones, which was presented in the previous chapter, and the prediction of the Lorenz attractor, which is a dynamic system known for its chaotic behaviour. Some comparisons are made with respect to some baseline techniques such as quantile regression and set-membership methods. In both examples, the proposed approach obtained better results. In the second half of the chapter, it was proposed a method that not only allows us to obtain prediction regions of multivariate systems or time-series but also reduces considerable the computational burden. This is due to the fact that the second method no longer needs the integration of the empirical conditioned density function. Also, a method to obtain ellipsoidal approximations of such regions was proposed. Finally, a prediction example of a kite system was proposed. Other approaches including GPs, IR, etc. were considered as baselines. It was shown that the proposed approach and the approximated ellipsoidal regions achieved the smallest size while fulfilling the desired probabilistic specifications. Part II Kriging-based identification 55 Chapter 4 State-space kriging for autonomous systems In part I, some methods to make predictions of time-series or dynamical systems based on dissimilarity functions were proposed. These included l-steps ahead forecastings, interval predictions, predictions regions, etc. From now on, the objective will be to obtain a dynamic model of a system given some past data of its outputs and inputs. In the case of autonomous systems only past outputs will be considered, whereas for non-autonomous systems past inputs will be considered as well. The state-space of this new model will be compounded of the vector of weights λappearing in the dissimilarity function presented in part I. Two different strategies are proposed, a linear time variant approach based on the weighting of the local data within the dissimilarity function and a kernel-based approach where the dissimilarity function is slightly modified to accommodate the kernel trick. Finally, some numerical examples are presented to show the effectiveness of the proposed approaches. 4.1 Dynamic kriging Consider an autonomous discrete nonlinear system xk+1 =h(xk) (4.1a) yk=g(xk),(4.1b) where kis the time instant, xk∈Rnxis the state of the system, yk∈Rnyis the output of the system, whereas h(·) and g(·) are unknown nonlinear functions such that h(·) : Rnx→Rnxand g(·) : Rnx→Rny. The objective of this section is to obtain a model of the outputs of (4.1). It is assumed that the only available data are the measurable outputs. The regressor zkis defined as a time delay embedding vector zk∈Rnzcontaining the 56 Chapter 4. State-space kriging for autonomous systems nppast outputs of the system. That is, zk= [y⊤ k, y⊤ k−1, ..., y⊤ k−np+1]⊤∈Rnz, where nz=npny. Also, note that z+ kis denoted as the successor of zk, i.e. z+ k= [y⊤ k+1, y⊤ k, ..., y⊤ k−np+2]⊤∈Rnz. From now on, assume that some historical data of the plant is stored in a database in the form of matrices: D=¯z1¯z2. . . ¯zN, D+=¯z+ 1¯z+ 2. . . ¯z+ N, where N > nzis the number of data points, ¯zrefers to a sample of zand matrix D+is the successor of D. The indexes of the columns of Dand D+do not refer to the sample time, but to the position in the matrix. Therefore, ¯zi+1 is not necessarily the successor sample of ¯zi. At sample time k, an estimation of the successor of zk, denoted as ˜zk+1, can be obtained by a linear combination of the columns of matrix D+using the vector of optimal weights λ∗ k∈RNobtained from an optimization problem similar to those shown in section 2, i.e. ˜zk+1 =D+λ∗ k. Considering definition 2.1 and assuming that γ= 0, the optimization problem from which λ∗ kis obtained can be posed as λ∗ k= arg min λk λ⊤ kH1λk(4.2a) s.t. D 1λk=zk 1,(4.2b) where H1∈RN×Nis a positive definite weighting matrix and 1a row vector with all its components equal to 1. Forcing the components of λkto sum one is equivalent to including a bias term in the estimation process. The simplest choice is to make H1equal to the identity matrix IN(see remark 2.2). This optimization problem can be rewritten as λ∗ k= arg min λk λ⊤ kH1λk(4.3a) s.t. Cλk=b , (4.3b) with C=D 1, b =zk 1. In order to guarantee that any point in Rnz+1 can be expressed as a linear combination of the columns of C, it is assumed that matrix Cis full row rank. This equality constrained quadratic problem has an analytic solution that can be obtained computing the Lagrangian and its derivative: L(λk, ν) = λ⊤ kH1λk+ν⊤(Cλk−b) 4.2. Linear state-space kriging 57 d dλL(λk, ν) = 2H1λk+C⊤ν , where νis the dual variable associated with the equality constraint. From the Karush-Kuhn-Tucker (KKT) conditions: 2H1λ∗ k+C⊤ν∗= 0 ,(4.4) which leads to λ∗ k=−H−1 1C⊤ν∗ 2. Pre-multiplying this equality by C, and taking into account that Cλ∗ k=b, the following is obtained ν∗=−2CH−1 1C⊤−1b, which applied to equation (4.4) yields λ∗ k=H−1 1C⊤CH−1 1C⊤−1zk 1. In order to predict zk+d, with d > 1, one could use this approach in a recursive way. That is, the i-th ahead prediction ˜zk+icould be used to compute λ∗ k+i=H−1 1C⊤CH−1 1C⊤−1˜zk+i 1, and thus obtaining ˜zk+i+1 =D+λ∗ k+i. In the next section, it is proposed a modification of this naive recursive method. The novel methodology relies on a timevarying state-space modelling of the optimal weighting vector parameter λ∗ k. 4.2 Linear state-space kriging Suppose that the prediction ˜zk+1 of zk+1 is obtained from ˜zk+1 =D+λ∗ k, where the sum of the components of λ∗ kis assumed to be equal to one. In order to model how the dynamics of the optimal vector of weights λ∗ k+1 depends on λ∗ k, a regularization term is added to optimization problem (4.2) that penalizes the difference between λ∗ k+1 and λ∗ k. In this way, vector λ∗ k+1 not only fulfills the required equality constraints, but also does not depart excessively from λ∗ k. This will reduce the sensitivity to noise of the identified dynamics. Thus, given ˜zk+1, λ∗ k+1 is obtained from λ∗ k+1 = arg min λk+1 (λk+1 −λ∗ k)⊤H2(λk+1 −λ∗ k) + λ⊤ k+1H1λk+1 s.t. D 1λk+1 =˜zk+1 1, 58 Chapter 4. State-space kriging for autonomous systems where H2∈RN×Nis chosen as the identity matrix multiplied by a certain scalar that could be selected by cross-validation [1, §16.5]. Because of the assumptions on λ∗ k, the previous optimization problem can be rewritten as λ∗ k+1 = arg min λk+1 (λk+1 −λ∗ k)⊤H2(λk+1 −λ∗ k) + λ⊤ k+1H1λk+1 s.t. D 1λk+1 =D+ 1λ∗ k. Thus, the problem becomes λ∗ k+1 = arg min λk+1 (λk+1 −λ∗ k)⊤H2(λk+1 −λ∗ k) + λ⊤ k+1H1λk+1 (4.5a) s.t. Cλk+1 =C+λ∗ k,(4.5b) with C=D 1, C+=D+ 1. Note that λ∗ k+1 is determined only by λ∗ kand matrices Cand C+. Optimization problem (4.5) can be rewritten as λ∗ k+1 = arg min λk+1 1 2λ⊤ k+1Hλk+1 +f⊤λk+1 s.t. Cλk+1 =b , with H= 2(H1+H2), f=−2H2λ∗ kand b=C+λ∗ k. Also, note that the constant term λ∗ k ⊤H2λ∗ kis removed because it does not affect the solution λ∗ k+1. The Lagrangian of this problem is given by L(λk+1, ν) = 1 2λ⊤ k+1Hλk+1 +f⊤λk+1 +ν⊤(Cλk+1 −b), where νis the dual variable associated with the equality constraint. The derivative of the Lagrangian is d dλk+1 L(λk+1, ν) = Hλk+1 +f+C⊤ν . In the optimum, the derivative fulfills the KKT conditions [139, §10.1.1]. That is, Hλ∗ k+1 +f+C⊤ν∗= 0 , and thus λ∗ k+1 =−H−1f−H−1C⊤ν∗.(4.7) Pre-multiplying both sides of last equality by Cyields b=Cλ∗ k+1 =−CH−1f−CH−1C⊤ν∗, 4.2. Linear state-space kriging 59 and thus ν∗= (CH−1C⊤)−1(−CH−1f−b). Substituting this into equation (4.7), the following expression for λ∗ k+1 is obtained λ∗ k+1 =H−1C⊤CH−1C⊤−1CH−1f+b−H−1f. Taking into account that f=−2H2λ∗ kand b=C+λ∗ k, the state-space equation of λ∗ kcan be written as λ∗ k+1 =Aλ∗ k, with A= 2H−1H2+H−1C⊤(CH−1C⊤)−1(C+−2CH−1H2). Note that this means that λkfollows linear dynamics. Thus, a new model for the outputs of system (4.1) has been obtained using historical data of these outputs. This new autonomous system allows us to compute the next values of λand z. Note that, as only the first term of zkis needed (i.e. the term corresponding to yk), the model can be posed as λ∗ k+1 =Aλ∗ k yk=Y λ∗ k, where Ydenotes a matrix compounded of only the first rows of D, that is, Yis a matrix containing the samples ¯yi. 4.2.1 Initial condition The initial vector of optimal weights λ∗ 0is computed by means of the optimization problem in (4.2) substituting kby 0, that is λ∗ 0= arg min λ0 λ⊤ 0H1λ0 s.t. D 1λ0=z0 1. 4.2.2 Local-data approach Note that the previous model is linear and time-invariant as the matrix Ais constant. However, it is possible to weight the points in the data set with respect to z, which would encourage the use of local data and thus provide better results when identifying nonlinear systems. This can be done by choosing H1appropriately. In remark 2.2, it was shown that it is possible to consider a vector ωwithin the dissimilarity function to weight different elements in the data set. Here, the squared Euclidean distance is chosen to measure the dissimilarity of a certain zto each point of the data set. That is, 60 Chapter 4. State-space kriging for autonomous systems ω(z) =    (z−¯z1)⊤(z−¯z1) . . . (z−¯zN)⊤(z−¯zN)   ,(4.8) and thus H1= diag(ω(z)) , where diag(·) denotes a diagonal matrix RN×RNwhose non-zero entries are the components of the input vector. Note that, because of this change, this vector ω(z) needs to be computed at each time instant k, leading to a different matrix Afor every k, and thus leading to LTV dynamics for λ, that is λ∗ k+1 =Akλ∗ k yk=Y λ∗ k. 4.3 Kernel-based state-space kriging Kernels functions are widely used in the machine learning field. For example, Support Vector Machines (SVM) are supervised learning models that classify linearly separable data. That is, given a cloud of data points (each point belonging to a certain class), the problem of classification is defined as finding an hyperplane that divide the data into two sets. However, when considering nonlinear relations within the data, it is not possible to obtain an hyperplane that separates the data. In order to be able to classify this kind of data, it is needed to project this data into a high-dimensional space where it may become linearly separable. Thanks to the so-called kernel trick (see chapter 3 in [150]), it is possible to operate in a high-dimensional feature space without computing the coordinates of such space. Instead, it is only needed to compute the inner product of the images of all pairs of the data samples, which bypasses the computation of the coordinates of the feature space that may be cumbersome or even impossible. Thus, instead of using local data as it was presented in the previous section, one could resort to the use of kernels in order to model the nonlinear dynamics of the system. Also, the use of kernels will allow us to compute the matrices of the system only once, unlike the local data approach where the matrices needed to be recomputed at each sample time. In order to include the kernels, the optimization problem is modified, becoming λ∗ k+1 = arg min λk+1 (λk+1 −λ∗ k)⊤H2(λk+1 −λ∗ k) + λ⊤ k+1H1λk+1 +     N X i=1 φ¯ziλk+1,(i)− N X i=1 φ¯z+ iλ∗ k,(i)     2 Σ−1 φ s.t. 1λk+1 = 1 , 4.3. Kernel-based state-space kriging 61 where φ(·) : Rnz→ H refers to a nonlinear operator that maps Rnzinto a probably high dimensional space H,φ¯ziand φ¯z+ idenote φ(¯zi) and φ(¯z+ i) respectively and Σφis a positive definite matrix of appropriate dimensions. In what follows, it is shown that we do not need a precise knowledge of φ(·) to compute λk+1 as it suffices to compute, for a given pair a∈Rnzand b∈Rnz, the product ⟨φa, φb⟩=φaΣ−1 φφb. Note that the previous linear hard constraint on zkhas been changed to a penalty term on a high dimensional space by means of the kernel trick as it is no longer a linear constraint. Assuming that some data sets of the evaluation of φ(·) over the time delay embeddings of Dand D+are available (which due to the kernel trick are not really necessary), it would be possible to denote them as φ¯z=φ¯z1φ¯z2. . . φ¯zN, φ¯z+=hφ¯z+ 1φ¯z+ 2. . . φ¯z+ Ni. Thus, the previous problem could be written in matrix form as λ∗ k+1 = arg min λk+1 (λk+1 −λ∗ k)⊤H2(λk+1 −λ∗ k) + λ⊤ k+1H1λk+1 +∥φ¯zλk+1 −φ¯z+λ∗ k∥2 Σ−1 φ s.t. 1λk+1 = 1 . Now, operating with the term ∥φ¯zλk+1 −φ¯z+λ∗ k∥2 Σ−1 φ, the following expression is obtained ∥φ¯zλk+1 −φ¯z+λ∗ k∥2 Σ−1 φ=λ⊤ k+1 φ⊤ ¯zΣ−1 φφ¯zλk+1 −2λ⊤ k+1 φ⊤ ¯zΣ−1 φφ¯z+λ∗ k +λ∗ k ⊤φ⊤ ¯z+Σ−1 φφ¯z+λ∗ k. Again, the constant term λ∗ k ⊤φ⊤ ¯z+Σ−1 φφ¯z+λ∗ kis discarded because it does not affect the values of λ∗ k+1, leading to the following optimization problem λ∗ k+1 = arg min λk+1 1 2λ⊤ k+1Hλk+1 +f⊤λk+1 (4.12a) s.t. 1λk+1 = 1 .(4.12b) with H= 2(H1+H2)+2φ⊤ ¯zΣ−1 φφ¯z, f =−2(H2+φ⊤ ¯zΣ−1 φφ¯z+)λ∗ k. Note that the kernel related terms can be computed because only cross products appear in the aforementioned equations. These terms would be computed as φ⊤ ¯zΣ−1 φφ¯z=     ⟨φ¯z1, φ¯z1⟩ ⟨φ¯z1, φ¯z2⟩. . . ⟨φ¯z1, φ¯zN⟩ ⟨φ¯z2, φ¯z1⟩ ⟨φ¯z2, φ¯z2⟩. . . ⟨φ¯z2, φ¯zN⟩ . . .. . .. . . ⟨φ¯zN, φ¯z1⟩ ⟨φ¯zN, φ¯z2⟩. . . ⟨φ¯zN, φ¯zN⟩      , 68 Chapter 4. State-space kriging for autonomous systems 0 5 10 15 -5 0 5 o K-SSK Real Measured Output 0 5 10 15 -5 0 5 p 0 5 10 15 Time (s) -5 0 5 l Figure 4.2: Forecasting the noisy R¨ossler attractor with kalman filtering. LD-SSK K-SSK GP NARX K-DMD RC 62.909 37.5912 2.591 14.872 0.036 0.868 Table 4.3: Average online computational time in milliseconds (R¨ossler attractor). 4.6 Conclusions This chapter presented the state-space kriging method for autonomous systems. The proposed technique allowed us to obtain a model of the system by means of past data of the process. This could be done by manipulating appropriately the dissimilarity function presented in part I and by obtaining the explicit solution of the resulting optimization problems. Also, it was shown that a kalman filter can be used to improve the performance of the proposed models. Finally, two numerical examples were proposed in order to show the effectiveness of the SSK approaches. Both the LD-SSK and the K-SSK attained better performance than other existing machine learning models. In the next chapter, the proposed methodology will be extended to tackle nonautonomous system, which will allow us to develop MPC controllers as well. 69 Chapter 5 State-space kriging for non-autonomous systems The previous chapter introduced the concept of state-space kriging for autonomous systems and, by means of numerical results, it was shown that the performance is, if not better, comparable to many other forecasting methods. However, in order to be able to use the proposed approaches in a control scheme, it is needed to introduce the input term in the SSK formulation. In this chapter, both methodologies presented before will be extended to tackle non-autonomous systems. It is also shown that the kernel-based SSK is more suitable for control due to the fact that the system matrices are computed only once, unlike the local data approach. 5.1 Non-autonomous linear SSK Consider a non-autonomous discrete nonlinear system xk+1 =h(xk, uk) (5.1a) yk=g(xk),(5.1b) where kis the time instant, xk∈Rnxis the state of the system, uk∈Rnuis the input of the system, yk∈Rnyis the output of the system and h(·) and g(·) are unknown nonlinear functions such that h(·) : Rnx×nu→Rnxand g(·) : Rnx→Rny. Unlike in chapter 4, it is assumed that the vector of weights λ∗ kfulfill the following constraints: zk=D+λ∗ k, 1 = 1λ∗ k. 70 Chapter 5. State-space kriging for non-autonomous systems Furthermore, it is imposed that λ∗ k+1 must be able to compute zkand ukas an affine combination of Dand U. Thus, it must satisfy the constraints   D 1 U λ∗ k+1 =  zk 1 uk . where U=¯u1¯u2. . . ¯uN∈ Rnu×Nis the data set of control actions, i.e. ¯ui are past samples of the input u. This set of constraints is similar to the constraints presented in the previous chapter but aimed to obtain the successor of λ∗ k+1, as a function of λ∗ kand uk, that is, to consider non-autonomous systems. Thus, the proposed optimization problem leading to the state equation is λ∗ k+1 = arg min λk+1 (λk+1 −λ∗ k)⊤H2(λk+1 −λ∗ k) + λ⊤ k+1H1λk+1 (5.3a) s.t. Cλk+1 =C+λ∗ k+0(nz+1)×nu Inuuk,(5.3b) where 0(nz+1)×nuis a matrix of zeros compounded of nz+ 1 rows and nucolumns C=D⊤1⊤U⊤⊤and C+=hD+⊤1⊤0⊤i⊤.Note that ukis the current value of the system input, not to be confused with ¯u, which are past values of ukstored in the data base U. This problem can be written in canonical form as λ∗ k+1 = arg min λk+1 1 2λ⊤ k+1Hλk+1 +f⊤λk+1 s.t. Cλk+1 =b , with H= 2(H1+H2), f=−2H2λ∗ kand b=C+λ∗ k+0nu×(nz+1) Inu⊤uk. Note that the constant term is discarded because it does not affect the value of λ∗ k+1. By means of the KKT conditions, this problem has the solution λ∗ k+1 =−H−1f+H−1C⊤CH−1C⊤−1CH−1f+b. Substituting fand b, the following state-space equation is obtained λ∗ k+1 =Aλ∗ k+Buk, yk=Y+λ∗ k. where Y+denotes a matrix compounded of only the first nyrows of D+and A= (2H−1H2+H−1C⊤(CH−1C⊤)−1C+−2CH−1H2), B=H−1C⊤(CH−1C)−10nu×(nz+1) Inu⊤. As the expression of the dynamics of λhas been obtained, the output ycan be easily obtained for any time instant kassuming that the inputs are given. 5.2. Non-autonomous kernel-based SSK 71 5.1.1 Initial condition It still remains to show how to obtain the initial value λ∗ 0. For that purpose, consider the following optimization problem λ∗ 0= arg min λ0 λ⊤ 0H1λ0 s.t. D+ 1λ0=z0 1. where z0is the initial value of the time delay embedding. 5.2 Non-autonomous kernel-based SSK Taking into account kernel functions, the optimization problem in (5.3) becomes the following λ∗ k+1 = arg min λk+1 (λk+1 −λ∗ k)⊤H2(λk+1 −λ∗ k) + λ⊤ k+1H1λk+1 +∥φ¯zλk+1 −φ¯z+λ∗ k∥2 Σ−1 φ(5.6a) s.t. U 1λk+1 =uk 1.(5.6b) As in the previous chapter, the hard constraint becomes a penalty term in a high dimensional space due to the kernel trick [150]. Now, operating with the term ∥φ¯zλk+1 −φ¯z+λ∗ k∥2 Σ−1 φ, the following expression is obtained ∥φ¯zλk+1 −φ¯z+λ∗ k∥2 Σ−1 φ=λ⊤ k+1φ⊤ ¯zΣ−1 φφ¯zλk+1 −2λ⊤ k+1φ⊤ ¯zΣ−1 φφ¯z+λ∗ k +λ∗ k ⊤φ⊤ ¯z+Σ−1 φφ¯z+λ∗ k, Discarding the constant term λ∗ k ⊤φ⊤ ¯z+Σ−1 φφ¯z+λ∗ k, we obtain the following optimization problem λ∗ k+1 = arg min λk+1 1 2λ⊤ k+1Hλk+1 +f⊤λk+1 s.t. Tλk+1 =0⊤ 1+Inu 0uk. with H= 2(H1+H2)+2φ⊤ ¯zΣ−1 φφ¯z, f =−2(H2+φ⊤ ¯zΣ−1 φφ¯z+)λ∗ k,T=U 1. Note that the kernel related terms can be computed because only cross products appear in the aforementioned equations. These terms would be computed as φ⊤ ¯zΣ−1 φφ¯z=     ⟨φ¯z1, φ¯z1⟩ ⟨φ¯z1, φ¯z2⟩. . . ⟨φ¯z1, φ¯zN⟩ ⟨φ¯z2, φ¯z1⟩ ⟨φ¯z2, φ¯z2⟩. . . ⟨φ¯z2, φ¯zN⟩ . . .. . .. . . ⟨φ¯zN, φ¯z1⟩ ⟨φ¯zN, φ¯z2⟩. . . ⟨φ¯zN, φ¯zN⟩      , 72 Chapter 5. State-space kriging for non-autonomous systems φ⊤ ¯zΣ−1 φφ¯z+=      ⟨φ¯z1, φ¯z+ 1⟩ ⟨φ¯z1, φ¯z+ 2⟩. . . ⟨φ¯z1, φ¯z+ N⟩ ⟨φ¯z2, φ¯z+ 1⟩ ⟨φ¯z2, φ¯z+ 2⟩. . . ⟨φ¯z2, φ¯z+ N⟩ . . .. . .. . . ⟨φ¯zN, φ¯z+ 1⟩ ⟨φ¯zN, φ¯z+ 2⟩. . . ⟨φ¯zN, φ¯z+ N⟩       . Applying the KKT conditions, the solution of the optimization problem (5.6) is obtained; λ∗ k+1 =Aλ∗ k+Buk+c(5.8a) yk=Y+λ∗ k,(5.8b) with A= 2H−1(IN−T⊤(TH−1T⊤)−1TH−1)(H2+φ⊤ ¯z Σ−1 φφ¯z+), B=H−1T⊤(TH−1T⊤)−1Inu 0, c=H−1T⊤(TH−1T⊤)−10⊤ 1. Remark 5.1. Note that adding the input to the optimization problem is compatible with both the Local-data approach and the kernel-based approach, leading to a family of different predictors. However, in the local-data approach, it is needed to compute the matrices Aand Bat each time instant k, becoming cumbersome when considering MPC control problems. For that reason, in the following, it is only considered the kernel-based SSK. Remark 5.2. Note that the kalman filter for the non-autonomous SSK is equivalent to the one proposed in the previous chapter. In this case, equation (4.14) becomes ˜ λ∗ k+1 =Aˆ λ∗ k+Buk+c . (5.9) 5.2.1 Initial condition As in the non-autonomous linear SSK, an initial value λ∗ 0is needed. This initial condition can be obtained from the following optimization problem λ∗ 0= arg min λ0 λ⊤ 0H1λ0+∥φz0−φ¯z+λ0∥2 Σ−1 φ s.t. 1⊤λ= 1 . Operating with ∥φz0−φ¯z+λ0∥2 Σ−1 φ, we obtain λ∗ 0= arg min λ0 1 2λ⊤ 0Hλ0+f⊤λ0 s.t. 1⊤λ0= 1 , 5.3. Application to MPC 73 with H= 2H1+ 2φ⊤ ¯z+Σ−1 φφ¯z+, f =−2φ⊤ ¯z+Σ−1 φφz0, φ⊤ ¯z+Σ−1 φφ¯z+=      ⟨φ¯z+ 1, φ¯z+ 1⟩ ⟨φ¯z+ 1, φ¯z+ 2⟩. . . ⟨φ¯z+ 1, φ¯z+ N⟩ ⟨φ¯z+ 2, φ¯z+ 1⟩ ⟨φ¯z+ 2, φ¯z+ 2⟩. . . ⟨φ¯z+ 2, φ¯z+ N⟩ . . .. . .. . . ⟨φ¯z+ N, φ¯z+ 1⟩ ⟨φ¯z+ N, φ¯z+ 2⟩. . . ⟨φ¯z+ N, φ¯z+ N⟩       , φ⊤ ¯z+Σ−1 φφz0=      ⟨φ¯z+ 1, φz0⟩ ⟨φ¯z+ 2, φz0⟩ . . . ⟨φ¯z+ N, φz0⟩       . 5.3 Application to MPC Assuming that we have an arbitrary nonlinear system in a discrete time setting like the system in equation (5.1), then, the objective of an MPC controller is to steer the state xto an equilibrium point (xs, us) that it is called the reference of the controller. The control actions applied to the system to achieve this objective ukare computed by means of an optimization problem. The obtained optimal inputs depend on both the actual state of the plant and the optimization criteria. This step cost function is designed to penalize the deviation of the state and the inputs to the reference at each time instant kfor a certain prediction horizon Np. In order to ensure stability for the controller, a terminal cost function is usually added to the previously defined step cost function. Denoting V(·) as the total cost function, uas the sequence of inputs {ui}k+Np−1 i=k, xas the sequence of predicted states {xi|k}Np i=0,xi|kas the prediction of xk+imade at instant kand considering that the system may have constraints in the inputs and in the states so that x∈ X and u∈ U, the optimal control inputs u∗are obtained by solving the following optimization problem u∗= arg min x,uV(x,u) (5.11a) s.t. system model (5.11b) xi|k∈ X ∀i= 0, . . . , Np(5.11c) uk+i∈ U ∀i= 0, . . . , Np−1.(5.11d) From this sequence of optimal inputs u∗, only the first component will be applied to the system, discarding the rest as it is usual in MPC due to the receding horizon strategy. In this section, a tracking MPC controller [124, 125] that uses the K-SSK model is presented in this section. The main difference of the tracking MPC with respect to 74 Chapter 5. State-space kriging for non-autonomous systems traditional MPC is that the reference to be tracked is considered as an additional decision variable in the optimization problem, that is, it becomes an artificial reference. In order to enforce that the artificial reference converges to the target reference, an additional cost is added to the MPC cost function, penalising the deviation between them. Among the many advantages of this formulation, the fact that the recursive feasibility is guaranteed for any change of the desired reference and a significantly larger domain of attraction for short prediction horizons are probably the most important ones. The cost function in the proposed controller is compounded of three ingredients: •A step cost function ls(·,·) to penalize tracking error. Tracking error with respect to a certain reference is tackled by means of a change of variables ˘y= ˜y−ys, ˘u=u−uswhere ˜yis a prediction of y,ysis the artificial output reference and usis the artificial input reference. Here, a quadratic cost is considered. This cost penalizes the distance to the artificial input and output reference by means of some weighting matrices of appropriate dimensions Qand R(i.e. Q∈Rny×nyand R∈Rnu×nu), i.e. ls(˘y, ˘u) = ˘y⊤Q˘y+ ˘u⊤R˘u. •An offset function penalizing the difference between the artificial output ys and the target desired reference, lo(ys, r)=(ys−r)⊤O(ys−r). Under some mild conditions [124], it is proven that the artificial reference converges to the true reference as time goes by. •A weighted terminal cost function. This function measures the closeness of the terminal state λ∗ Np|kto the artificial steady state λ∗ s. Weighting appropriately the terminal cost allows us to omit the terminal equality constraint [125], simplifying the design of the controller. lt(λ∗ Np|k, λ∗ s) = γλ∗ Np|k−λ∗ s⊤Pλ∗ Np|k−λ∗ s, where γ≥1. Denoting y∈RnyNpand u∈Rnu(Np+1) as y= [y⊤ 0|k, . . . , y⊤ Np−1|k]⊤,[u⊤ k, . . . , u⊤ k+Np−1, u⊤ s]⊤, it is possible to define a total cost function VNp(y,u, r, λ∗ k) as the sum of the aforementioned three functions for a finite prediction horizon Np VNp(y,u, r, λ∗ k) = Np−1 X i=0 ls(˘yi|k,˘uk+i) + lo(ys(us), r) + lt(λ∗ Np|k(λ∗ k,u), λ∗ s(us)). 5.3. Application to MPC 75 Thus, an MPC control problem with inequality constraints in the outputs and box constraints in the inputs is presented here u∗= arg min y,u,λNp|k VNp(y,u, r, λ∗ k) (5.12a) s.t. λ∗ i+1|k=Aλ∗ i|k+Buk+i+c∀i= 0, . . . , Np−1 (5.12b) yi|k=Y+λ∗ i|k∀i= 0, . . . , Np(5.12c) λ∗ s=Aλ∗ s+Bus+c(5.12d) ψ yi|k≤δ∀i= 0, . . . , Np(5.12e) umin ≤uk+i≤umax,∀i= 0, . . . , Np−1 (5.12f) which is a parametric quadratic optimization problem whose parameters are r and λ∗ k. Appendix A shows how to pose this problem in canonical form. From this optimization problem, a sequence of optimal control actions uis obtained. However, due to the receding horizon scheme characteristic of any MPC controller, only the first component is applied to the system, computing a whole new sequence at next time instant. 5.3.1 Nominal stability analysis For the nominal stability analysis, it is assumed that there are no mismatches between the dynamics of the real system and the obtained prediction model. Furthermore, the state vector λ∗ kis assumed to be known. This implies that the kalman filter is not necessary, and thus it is not taken into account. First, note that the proposed model is not strictly linear with respect to the weights λk. However, taking into account that λ∗ k+1 =Aλ∗ k+Buk+c λ∗ s=Aλ∗ s+Bus+c, and subtracting the equations, it is clear that (λ∗ k+1 −λ∗ s) = A(λ∗ k−λ∗ s) + B(uk−us). Making ˘ λ∗ k=λ∗ k−λ∗ s,˘uk=uk−us, it is easy to see that the dynamics of ˘ λ∗ kare linear, that is ˘ λ∗ k+1 =A˘ λ∗ k+B˘uk. Now, some additional assumptions are made in order to prove the nominal stability of the controller: 76 Chapter 5. State-space kriging for non-autonomous systems Assumption 5.3. The system described in (5.8) is observable and the pair (A, B) is stabilizable. Assumption 5.4. Q,Rand Oare positive definite matrices. Assumption 5.5. It exists a stabilizing matrix Kso that the matrix (A+BK) is Hurwitz. Assumption 5.6. It exists a matrix Pso that the Lyapunov equation (A+BK)⊤P(A+BK)−P=−(Q+K⊤RK) (5.13) is fulfilled. Assumption 5.7. It exists at least one feasible equilibrium point satisfying the constraints of the controller. The previous assumptions makes possible to establish the following closed-loop stability lemma. Lemma 5.8 (Closed-loop stability).The system (5.8) controlled with the proposed MPC controller (5.12) converges asymptotically to the desired reference. In the case where the desired reference is unreachable, the system converges to the feasible reference which minimizes the offset cost while maintaining the asymptotic stability. Proof. Assuming that the assumptions are satisfied, the proof follows from [124] and [125]. In [124], it is proven the asymptotic stability of the tracking controller for linear systems using a terminal cost term and a terminal constraint. Then, in [125], it is proven that the stability can be maintained even if there is no terminal constraint. In that case, increasing the weighting of the terminal cost makes the domain of attraction to grow. ■ 5.3.2 Robust stability analysis In this section, it is considered that modeling errors may appear, due to the fact that the model of the system is not perfect. For this purpose, the Robustly asymptotically stability (RAS) notion showcased in [106] is chosen. Denote eas measurement errors and das additive disturbances. Also, denote eand das sequences of Npelements of eand d. Then, the definition of RAS is the following. 5.4. Examples 77 Definition 5.9. The origin of the closed loop nonlinear system xk+1 = h(xk, κ(xk)) is considered to be RAS in the interior of a set Fwith respect to both measurement errors and additive disturbances if it is possible to find a KL function β[152] and for each ϵ > 0 and compact set C ⊂ F there exists δ > 0 such that 1. max(d)≤δ, max(e)≤δ. 2. Any admissible trajectory is bounded by β(|x|, k) + ϵ. Lemma 5.10 (Robust asymptotically stability).The system (5.8) joined together with the kalman filter observer and controlled with the proposed MPC controller (5.12) is RAS with respect to measurement noise and additive disturbances. Proof. Proposition 8 in [106] establishes some sufficient conditions to guarantee RAS. Actually, it is only needed to prove the continuity of the Lyapunov function. Here, in the proposed controller, taking into account that the prediction model is linear, then, the parametric quadratic optimization problem is convex and the aforementioned condition holds. Thus, the system is RAS with respect to measurement noise and additive disturbances. ■ 5.4 Examples This section presents two application examples of the proposed controller. First, the application to a simulated single-input single-output system, a continuously stirred tank reactor, is presented, followed by the application to a laboratory temperature control equipment. This is a multivariable system with two inputs and two outputs. Some considerations are taken into account for both examples in order to guarantee that the assumptions are fulfilled: 1. The models obtained for both systems are open-loop stables, i.e. all the eigenvalues in Aare less than the unity. Thus, their modes are stable and hence the system is stabilizable (assumption 5.3). 2. Matrices Q,Rand Oare chosen as the identity matrix of appropriate dimensions times some positive constant, thus they are positive definite (assumption 5.4). 3. Kand Pare computed using the Matlab “dlqr” function. This returns a stabilizing gain Kand the matrix Pthat satisfies the Lyapunov equation (5.13) (assumptions 5.5 and 5.6). 4. Taking into account the previous statements and that the following examples only consider box constraints in the inputs, it is easy to see that assumption 5.7 is also fulfilled. 84 Chapter 5. State-space kriging for non-autonomous systems Once reached this point, it is possible to write the output tracking error term in matrix form. Now, making L=   Inu0 0 −Inu . . ..... . .. . . 0 0 Inu−Inu   ,Q=   Q0 0 0...0 0 0 Q   ,R=   R0 0 0...0 0 0 R   , we obtain Np−1 X i=0 ls(˘yi|k,˘uk+i)=(y−ys)⊤Q(y−ys) + u⊤L⊤RLu = (Mu +n−WYu+vY)⊤Q(Mu +n−WYu+vY) + u⊤L⊤RLu. Also, lo(ys, r) can be written as lo(ys, r) = Y+Wu+Y+v−r⊤OY+Wu+Y+v−r. It only remains to obtain a matrix expression for lt(λ∗ Np|k, λ∗ s). As it was already shown how to obtain λ∗ s, it is only needed to show how to obtain λ∗ Np|k. Using equation (5.15), it is clear that λ∗ Np|k=Gu+h, where G=ANp−1B ANp−2B··· B0, h=ANp−1λ∗ k+ Np−1 X i=0 ANp−1−ic, and thus lt(λ∗ Np|k, λ∗ s) = γ(Gu+h−Wu−v)⊤P(Gu+h−Wu−v). It is easy to see that the optimization problem in (5.12) becomes min u 1 2u⊤Hu+f⊤u+ cnt s.t. ΨMu ≤∆−Ψn u≤u≤u. where H= 2((M−WY)⊤Q(M−WY) + L⊤RL +W⊤Y+⊤OY +W +γ((G−W)⊤P(G−W))), f⊤= 2(n−vY)⊤Q(M−WY) + 2(Y+v−r)⊤OY +W+ 2γ(h−v)⊤P(G−W), cnt = (n−vY)⊤Q(n−vY)+(Y+v−r)⊤O(Y+v−r) + γ(h−v)⊤P(h−v). 5.5. Conclusions 85 Thus, the control action is computed as u∗= arg min u 1 2u⊤Hu+f⊤u s.t. ΨMu ≤∆−Ψn u≤u≤u, where only the first component of the minimizer sequence u∗is applied following a receding horizon scheme as usual in any predictive controller. 86 Chapter 5. State-space kriging for non-autonomous systems Part III Probabilistically-certified data center management 89 Chapter 6 Bounds on the constraint violation level 6.1 Introduction Bounds on the constraint violation levels are useful in constrained control problems as they provide a measure on how likely are the constraints to be violated in practice. Being MPC the most popular control technique using constraints, it will be chosen for the example in this chapter. Consider a nonlinear system which is also affected by noises and disturbances, that is xk+1 =h(xk, uk, wk) yk=g(xk, vk), where kis the time instant, xk∈Rnxis the state of the system, uk∈Rnuis the input of the system, yk∈Rnyis the output of the system, wkare disturbances in the state, vkis measurement noise and h(·) and g(·) are nonlinear functions such that h(·) : Rnx×nu→Rnxand g(·) : Rnx→Rny. In this case, the control action computed using a deterministic MPC law may not be optimal for every possible realization of wkand vk. Also, it is not possible to guarantee that the constraints will be satisfied. For that reason, it may be needed to resort to robust MPC schemes, which are known to be overly conservative in most of the cases as it was showcased in the introduction. Another possibility would be to rely on strategies based on randomized settings such as stochastic MPC or chance-constrained MPC, which may reduce this conservatism. However, this may become an intractable problem sometimes. Instead, in this chapter, two methodologies to bound the constraint violation rate of a finite family of controllers by means of offline simulations of the closed-loop 90 Chapter 6. Bounds on the constraint violation level system are presented. This finite family of controllers can be compounded of any controller, not necessarily MPC controllers. The obtained bounds will allow us to compare the proposed controllers in order to obtain the best one according to a pre-specified criterion that could weight performance and constraint violation. Also, it does not increase the computation burden of the control problem online. 6.2 Main results Assume that a certain controller is parameterized by means of a decision vector θ∈Θ, where Θ is a set compounded of all the possible values that the vector θcan take. The disturbance wfollows a certain probability distribution whose sample space is W. No assumptions are made with respect to the size or shape of W. From this, a function q: Θ ×W → {0,1}is defined. This will be helpful to formulate the problems of this section in a general setting. Thus, the function q(·) is as follows q(θ, w) = (0 if θfulfills some design specifications for w 1 otherwise. This specification could be, for example, the satisfaction of a certain constraint subject to the realization w. Then, given θ, the probability of violation of the aforementioned constraint for any w∈ W can be written as ψ= Prob {q(θ, w) = 1 |θ}. As it is almost impossible to compute or measure this quantity exactly, we rely on the empirical mean to approximate this value. Consider wias the i-th realization of the disturbance w, and assuming a number of Nsamples, then ψ=1 N N X i=1 q(θ, wi), where denotes the empirical mean. Thus, ψwould be the empirical mean of ψfor any given experiment. Note that as ψis also a random variable, any set of experiments, that is, Ndifferent realizations of q(θ, wi) lead to different values for ψ. The objective of this chapter is to obtain tight bounds for this empirical violation rate. Definition 6.1. Failure. Assuming that E(ψ)> ρ + ∆ρ, where E(·) is the expectation operator, ρ∈[0,1] and ∆ρ > 0. A failure is considered to happen when an empirical violation rate of the controller ψfulfills ψ≤ρ. 6.2. Main results 91 That is, a failure happens when the empirical violation rate leads to misleading results, i.e., when given a empirical violation rate ψ≤ρ, the real violation rate ψ increases a quantity larger than ∆ρ. The following theorems prove that, with a properly chosen ∆ρ, failures occur with a probability lower than a small confident parameter δ. This definition of failure is slightly different than the one presented, for example, in [158]. Then, the following theorems summarize the results of this chapter. 6.2.1 A first bound on constraint violation rate Theorem 6.2 (Bound on the empirical constraint violation level).Assuming that E(ψ)> ρ + ∆ρ , with ρ∈[0,1] and ∆ρsatisfying ∆ρ≥1 Nlog 1 δ+ 2rρ Nlog 1 δ.(6.1) Then, Prob ψ≤ρ< δ , which implies that the real constraint violation level is bounded by ρ+ ∆ρwith probability greater than 1 −δ. Proof. The objective is to upper bound Prob ψ≤ρunder the assumption that E(ψ)> ρ + ∆ρ. Since the number of empirical violations grows with E (ψ), it is clear that this probability is lower bounded by the probability when E (ψ) = ρ+ ∆ρ. That is, Prob ψ≤ρ|E (ψ)> ρ + ∆ρ<Prob ψ≤ρ|E (ψ) = ρ+ ∆ρ. Once reached this point, it is easy to notice that the probability of a certain constraint to be violated can be interpreted as a Bernoulli random variable. Therefore, the probability of having less than a certain number of violations for some number of trials can be expressed as the binomial tail. On the other hand, consider a random variable Xfollowing a binomial random distribution B∼(N, p) where Nis the number of trials and pis the probability of success. Then, given a constant x, it is obtained from Chernoff’s bound [128, 159] that, if x≤p, then Prob X N≤x|EX N=p≤e−Nφ(ρ,ρ+∆ρ), where φ(x, p) = xlog x p+ (1 −x) log 1−x 1−p. 92 Chapter 6. Bounds on the constraint violation level The convergence rate provided by this bound is known to be tight from Cram´er’s theorem of large deviations when N→ ∞ ([160], chapter 23 in [161]). It is easy to see that the above Chernoff’s bound can be applied to our problem making X N=ψand ∆ρ > 0, obtaining Prob ψ≤ρ|E (ψ) = ρ+ ∆ρ≤e−Nφ(ρ,ρ+∆ρ). Thus, this upper-bound can be designed to be lower than the confidence parameter δby imposing e−Nφ(ρ,ρ+∆ρ)≤δ . Taking logarithms in both sides and rearranging the terms φ(ρ, ρ + ∆ρ)≥1 Nlog1 δ.(6.2) As ∆ρis embedded within the function φ(ρ, ρ + ∆ρ), the bound presented by Okamoto (Lemma 2 in [159]) is applied to obtain a closed expression for ∆ρ, φ(ρ, ρ + ∆ρ)≥pρ+ ∆ρ−√ρ2. Thus, the following sufficient condition for equation (6.2) is obtained pρ+ ∆ρ−√ρ2≥1 Nlog1 δ. For ∆ρ≥0, this is equivalent to pρ+ ∆ρ≥√ρ+r1 Nlog1 δ, and thus ρ+ ∆ρ≥ √ρ+r1 Nlog1 δ!2 . This is easily rewritten as ∆ρ≥1 Nlog 1 δ+ 2rρ Nlog 1 δ. This completes the proof. ■ 6.2.2 A different bound A different bound will be given in the following. This new bound is not guaranteed to be tighter than the previous one but it can yield better results in some cases. 6.2. Main results 93 Theorem 6.3. Suppose that E(ψ)> ρ + ∆ρ , with ρ∈[0,1] and ∆ρ≥log 1 δ+⌊ρN⌋log a N1−1 a−ρ, ∀a≥1.(6.3) Then, Prob ψ≤ρ< δ , As in theorem (6.2), this implies that the real constraint violation level is bounded by ρ+ ∆ρwith probability greater than 1 −δ. Proof. Assuming that E(ψ)> ρ + ∆ρ , one could write this probability as Prob ψ≤ρ|E (ψ)> ρ + ∆ρ. Again, as the number of empirical violations grows with E (ψ), it is clear that this probability is lower bounded by the probability when E (ψ) = ρ+ ∆ρ, i.e. Prob ψ≤ρ|E (ψ)> ρ + ∆ρ<Prob ψ≤ρ|E (ψ) = ρ+ ∆ρ. Similarly, as in theorem 6.1, it is easy to notice that this probability can be seen as a Bernoulli random variable and thus the probability of having less than a certain number of violations for some number of trials can be expressed as the binomial tail B(N, η, m) = m X y=0 N yηy(1 −η)N−y, where B(N, η, m) represents the mass probability function of the Binomial distribution. From the definition of ψ, it is easy to see that ⌊ψN⌋=ψN whereas ⌊ρN⌋ ≤ ρN. Then, Prob ψ≤ρ|E (ψ) = ρ+ ∆ρ= Prob ⌊ψN⌋≤⌊ρN⌋|E (ψ) = ρ+ ∆ρ. Now, it is possible to find the equivalence between this probability and the binomial tail as Prob ⌊ψN⌋≤⌊ρN⌋|E (ψ) = ρ+ ∆ρ= B(N, ρ + ∆ρ, ⌊ρN⌋).(6.4) Lemma 1 in [61] shows that the binomial tail can be bounded by the following expression B(N, η, m)≤amη a+ 1 −ηN,∀a≥1. 100 Chapter 7. Energy-efficient management of data centers example, the paper [162] proposes a control architecture where both thermal and tasks management are taken into account. The data center is posed as a linear first-order continuous system in its thermal component and in the computational part. Task arrival rates are considered deterministic, thus tasks arrive at fixed intervals. A Quality of Service (QoS) constraint is used to establish stability limits to the system. On the other hand, in [171], several control policies assuming that the maximum temperature of the servers are soft constraints are proposed. Other techniques available are based on workload distribution with some thermal aware criteria [172, 173]. Also, a solution to the problem under the assumption that the re-circulation of hot air is at constant temperature can be found in [174]. Model Predictive Control [86] has also been applied to this problem. In [175], a scheduling methodology based on MPC and electricity prices is proposed in order to reduce the economic impact. On the other hand, [176] unifies the management with the cooling scheme, assuming implicitly that servers have no limits to be overclocked, which is rather unrealistic. Due to the complexity of the problem, in [177] it is provided an approximation algorithm focused on providing fast evaluations of the complex constraints that have to be taken into account. However, the models appearing in the literature, under different assumptions, usually evade the fact that the system is actually a queue model, simplifying drastically the problem at the expense of a less realistic modeling. This chapter proposes an MPC framework for the optimization of cold aisle data centers. The optimization objective is to guarantee a certain QoS to the users and keep a suitable temperature of the servers while consuming the least amount of energy possible. The data center is modelled as a queue system where the arriving tasks can be computed by multiple servers at the same time. Although it moves away from the usual M/M/c queue and makes the treatment more complicated, it allows more generality within the model, where it is possible to take into account some specific kind of data centers like render farms [178]. Unlike other works available in the literature, the full queue model is used to predict the evolution of the system through the horizon, thus achieving a more realistic modeling. In order to tackle random arrival rates and workloads, two different strategies are proposed. First, it is considered that the arrival rates and workloads are deterministic, that is, it is assumed that they take always the value of the expectation of their respective random variables. On the other hand, a scenario-based approach where a certain number of realizations of the disturbances are drawn is proposed (see chapter 1). The QoS is managed by imposing a hard constraint on the number of powered servers in the optimization problem, thus guaranteeing that tasks are executed within the pre-specified time limits. Although there exist very optimized mixed-integer solvers nowadays such as Gurobi [179], they are not suitable for the related optimization problem. For that reason, particle algorithms were chosen to deal with the optimization problem [180]. These algorithms are highly parallelizable and can obtain important speed ups with a parallel implementation [181]. 7.2. Data center description 101 Finally, the bounds obtained in chapter 6 will be used to quantify the constraint violation rate of the proposed approaches. This makes possible to obtain a different characterization of the constraint violation rate for any proposed controller, which could be helpful to tune some parameters of the controllers like the control horizon or the number of scenarios to be considered within the optimization problem. 7.2 Data center description Data centers are compounded of thermally-isolated units in which server racks are allocated typically following a Cold Aisle (CA) structure as it can be seen in figure 7.1. Here, the air flow is separated into two different flows, being the first a cold one which reaches every server that it is blown from below the floor by means of suitable built-in fans. For this purpose, a Computer Room Air Conditioning (CRAC) unit is responsible of generating this cold flow. Then, the cold air travels through the servers and getting heated throughout this process and, finally, this hot air returns to the CRAC unit through the ceiling. CRAC CRAC Ceiling RACKRACK RACK CA Unit RACKRACK RACK Floor Hot air Cold air Figure 7.1: Scheme of a cold aisle data center structure. The CRAC unit is essential within the data center system because the server temperatures should be kept below certain security levels in order to ensure the servers reliability. As almost all power consumed by servers is dissipated as heat, CRAC operation is very costly because of the high number of servers that typical data center has [166]. Thus, any measure aimed to achieve an efficient management of the data center will have a great impact not only in operating costs and environmental impact, but also in the QoS provided to the clients. In the next sections, a discretized model of the data center dynamics with an integration time step of tsis presented. The discrete time unit will be denoted as k. 102 Chapter 7. Energy-efficient management of data centers 7.2.1 Tasks model In this section, the queue model of the data center operation is presented [171]. It is assumed that the data center has Mavailable servers and m, the number of booted servers, it is considered as an input of the system. In this way, the number of servers working in a certain time instant can be tuned according to the needs of the users (workload) and taking into account efficiency and QoS constraints. For the whole data center, there is a queue where the tasks wait until they can be processed. That is, until there is an available server. As it is accepted in the literature [171, 176], the time between arrivals is assumed to follow an exponential distribution with mean kaand a probability density function f(k) = 1 ka e−k ka.(7.1) Also, let Lkbe the request rate at instant k. That is, the number of tasks arriving to the data center at instant k. Because the time between arrivals was previously assumed to follow an exponential distribution, Lkcan be modeled as a Poisson random variable with mean 1 kaand probability mass function g(n) = 1 n!1 kan e−1 ka.(7.2) The number of tasks arriving in the interval [k, k −ku] is denoted as Lku,k for an interval of length ku>0, ku∈Z. Also, the probability mass function is gku(n) = 1 n!ku kan e−ku ka,(7.3) with mean ku ka. As it holds in practice, it is considered that tasks have a different workload. Then, for a certain task, the computational time required to complete it in a single server is assumed to follow an exponential distribution like (7.1) with mean 1 µ. Thus, Wkcorresponds to the average workload of the Lktasks that arrived at instant k. Similarly, Wku,k corresponds to the average workload of the Lku,k tasks that arrived between kand k−ku. Also, the parameters kaand µdefine the minimum number of servers that should be turned on within the data center so that the queue does not grow infinitely. Then, Mmust satisfy the condition M≥1 kaµ.(7.4) 7.2.2 Server model For a certain server i, the state xk,(i)at instant kis considered to be compounded of: 7.2. Data center description 103 1. The number of tasks currently running in the data center (αk). 2. The number of tasks in the queue (βk). 3. The temperature of the cold air (Tc,k). 4. The temperature of the server (Tk,(i)). 5. The time instant in which the server iis turned on (si). Obtaining the state of the whole data center is easy as it suffices to include the remaining Tk,(i)and s(i). Therefore, x⊤ kis obtained as x⊤ k=αk, βk, Tc,k, Tk,(1) . . . Tk,(M), s(1) . . . s(M). On the other hand, servers switch between four possible working conditions: •Off. The server does not draw any power. •Booting. This working condition appears due to the fixed delay kon that it takes to turn on a server from off to working or idle. While in this transition time, the server draws power at the same rate as the idle condition. •Working. The server is “on” and it is processing a certain task. For simplicity reasons, it is assumed that the power is drawn at a constant rate. In practice, this is equivalent to assume that the CPU frequency is constant. •Idle. The server is “on” but it is not processing any tasks, drawing less power than in the working state. However, the server is consuming energy for doing nothing. As usual in the literature and for simplicity reasons, the transition from “on” to “off” is assumed to be instantaneous. Also, the transition from working to idle and vice-versa is considered negligible. Note that because of the transition time kon and the power drawn throughout the booting process, it may be advantageous to keep a certain number of servers in idle state in the case when they are expected to be needed in a near future. Furthermore, the transition from “off” to “on” is always considered to be immediately available (although it takes a kon time to be completed), but the reverse, that is, from “on” to “off”, is done in a deferred way, because the server must finish the remaining tasks. Taken into account the previous working conditions, let uk,(i)be an input that indicates if the server iis switched on (uk,(i)= 1) or off (uk,(i)= 0) at a certain time instant k. Then, the power consumption of a server iis defined by the following conditions pk,(i)(xk,(i), uk,(i)) =            0 if off (uk,(i)= 0) a2if booting (uk,(i)= 1, αk≥0, k −si< kon) a2if idle (uk,(i)= 1, αk= 0, k −si≥kon) a1+a2if working (uk,(i)= 1, αk>0, k −si≥kon) 104 Chapter 7. Energy-efficient management of data centers where a1is the marginal consumption and a2the minimum consumption. The management approach presented in this thesis does not deal with the task scheduling of the data center [182]. For that reason, it is assumed that the distribution of the tasks among the servers follows some known rules. It is also assumed that a task can be split among multiple servers (up to M). However, a single server can work only for the completion of a single task. Because of the previous assumption, the data center will not work as a M/M/c queue in which each task is scheduled to be executed in a single server, complicating the mathematical modelling. However, it will result in a more general data center model. These server assignment policies imply the existence of a pool of running tasks whose length is equal to the number of servers in “on” condition (i.e., “working” or “idle”). Once a task is ready to be processed, that is, it is at the front of the queue and the pool has at least one empty slot, it is assigned to at least one server and it never leaves the pool until its completion. Thus, the remaining workload of each task is strictly decreasing. Denoting mkas the number of servers turned on at instant k(i.e. working or idle), the processing of a task is done in the following way. Assume that only a certain task is within a pool of mkservers (αk= 1). The workload of the task is the number of “work packages” that need to be processed to complete the task. At instant k,mkservers are assigned to this task (because it is the only task in the pool). Thus, every server will compute a “work package” resulting in mk “work packages” executed. If there are no “work packages” left for the instant k+ 1, the task is completed and ejected from the pool. If the pool is full at a certain instant k(i.e. αk=mk), the assignation is trivial because every task can only be assigned one server. For the case where 1 < αk< mk, the tasks will be assigned to one server and the task with largest remaining time will be assigned to mk−αkservers. This makes sense because it is easy to see that it will achieve better QoS. 7.2.3 Thermal model The thermal model of the servers is derived from the thermal balance equations Kt dT(i)(t) dt =cpqa(t)(Tc(t)−T(i)(t)) + p(i)(t),(7.5) where Tcand qaare the temperature and flow of the cold air provided by the CRAC unit, Ktthe server thermal capacity, T(i)and p(i)are the temperature and power consumption of server irespectively and cpthe air heat capacity. Considering a discretization scheme with an integration step ts, equation (7.5) turns into Tk+1,(i)=Tk,(i)+ts Ktcpqa,k Tc,k −Tk,(i)+pk,(i).(7.6) 7.2. Data center description 105 A hard constraint is used to guarantee the server reliability, keeping the temperatures under certain safety levels. That is T(i)≤80◦C∀i∈M . From the variables affecting T(i)in (7.6), the flow qa,k is assumed to be constant and only p(i)and Tccan be considered manipulable, the first one through the state of the server and the second one is assumed to be regulated by a set point Tr. That is, Tcfollows Trwith a first order closed loop dynamics with a time constant τand unity gain τdTc(t) dt =Tr(t)−Tc(t).(7.7) Similarly, equation (7.7) turns into Tc,k+1 =Tc,k +ts τ(Tr,k −Tc,k).(7.8) On the other hand, hard constraints in the input Trare considered 15◦C≤Tr≤25◦C. Also, the coefficient of perfomance (CoP) of the CRAC unit will change depending on the cold air temperature (Tc). The CoP represents how expensive is the cooling process of the air flow until a certain temperature. Usually, more power consumption is required to reach lower temperatures. As the CoP increases, the cost will decrease. It can be calculated from the following equation CoP(Tc,k) = 0.0068 T2 c,k + 0.0008 Tc,k + 0.458 ,(7.9) which is widely adopted in the literature [173]. Thus, let M P i=1 pk,(i)(xk,(i), uk,(i)) be the server power consumption at time instant k. The power consumption at the CRAC unit can be computed as M P i=1 pk,(i)(xk,(i), uk,(i)) CoP(Tc,k). Thus, the total power consumption corresponds to the power drawn by the servers added to the CRAC power consumption, leading to M X i=1 1 + 1 CoP(Tc,k)pk,(i)(xk,(i), uk,(i)). 106 Chapter 7. Energy-efficient management of data centers 7.2.4 Quality of service The QoS of a task is defined as the time required to finish it since its arrival until its completion. It includes both the waiting time within queue and the execution time once it is in the pool. Guaranteeing that this time will be lower than an agreed one is a necessary operating condition. In the case where the tasks are assigned to just one server, like in an M/M/c queue (or Erlang-C model [183]), the mean service time follows the expression tc,k =1 mk Wk−Lk .(7.10) This measurement could be used in practice with estimations of Wkand Lk, denoted as ˜ Wkand ˜ Lkrespectively. However, in the proposed approach, as a certain task can be executed concurrently in many servers, the aforementioned measure is used to provide an upper bound and thus, in this case, it can only be lower or equal tsv,k ≤1 mk ˜ Wk−˜ Lk , where ˜ Wkis the estimated value of Wkat instant kfor the next ˜ Lkrequests and ˜ Lkis the estimation of Lkat instant k. It is considered that the QoS is satisfied at instant kif the mean service time is not larger than a specified value D. In other words, the QoS constraint is satisfied if 1 mk ˜ Wk−˜ Lk≤D, which implies that mk≥˜ Wk1 D+˜ Lk. The term on the right represents the minimum number of servers to fulfill the QoS constraint with a certain ˜ L,˜ Wand D. This can be written as mD,k =˜ Wk1 D+˜ Lk. 7.3 Management approach In this chapter, an optimal management policy inspired on predictive control strategies is proposed. This controller decides the number of servers that should be on or off and the set point temperature of the CRAC unit. The optimization objective used to decide the optimal values of the inputs will be the minimization of the energy consumption and the control effort of the inputs (i.e. server switching and temperature reference changes), subject to thermal and QoS constraints. 7.3. Management approach 107 First of all, it is assumed that the dynamics of the queue, task arrivals, etc. run much faster than the process of switching on a server (i.e. kon ≫1). Having such a large delay in the control actions, the control decisions have to be made in a superior time scale and separate them over time. This leads us to a scheme where the predictive controller is not executed at every instant kbut it is executed at every instant kmtswhere kmis the number of instants between the controller execution. Thus, the sample time of the controller is kmts. In order to avoid switching off a server that is still booting (that is, never reached the “on” state and did nothing but consuming power), a sample time larger than the switching on time is chosen, i.e. km≥kon. Based on the model, the expected evolution of the data center can be estimated and then, it is possible to associate a predicted cost to a sequence of candidate future control inputs. For clarity purposes, denote ℓj=k+jkm. This will work as the time scale of the MPC controller. In this chapter, the following predicted cost function to measure the expected performance of the data center will be used: V(ℓ0,x,u,Tr) = Np X j=0 M X i=0 1 + 1 CoP(Tc,ℓj)pℓj,(i)(xℓj,(i), uℓj,(i)) +κu Nc X j=0|∆uℓj|+κTr Nc X j=0|∆Tr,ℓj|,(7.11) where Npand Ncare the prediction and control horizons, xis the sequence of xℓjover the prediction horizon, uand Trare the sequences of “on”-“off” control actions of all servers and temperature set points through the control horizon respectively, ∆uℓjis the total number of commutations to either “on” or “off” at instant ℓj, ∆Tr,ℓjis the increment in Tr,ℓj,κuis a term weighting ∆uℓjand κTris a term weighting ∆Tr,ℓj. Also, the control actions further from the control horizon are considered to remain constant. In order to derive the proposed controller, it is necessary to determine a prediction model such that for a given state at time ℓj,xℓjand for given control actions uℓjand Tr,ℓj(that will remain constant throughout the sampling time km), the state of the data center predicted at the next sampling time ˜xℓj+1 is calculated depending on the estimation of the number of tasks and their workload ˜ Lkm,ℓj+1 and ˜ Wkm,ℓj+1 . This prediction model can be posed as: ˜xℓj+1 =h(xℓj, uℓj, Tr,ℓj,˜ Lkm,ℓj+1 ,˜ Wkm,ℓj+1), being h(·) the function that compute the following state given the previous one, the inputs and the realisations of Lkmand Wkm. Note that the function h(·) must compute all events happening thorough the interval kmfor every instant k in order to be able to return the state at the following sample time. In practice, this function is evaluated by means of an open-loop simulation of the complete queue model of the data center. 108 Chapter 7. Energy-efficient management of data centers The optimal predicted number of servers and set point temperatures for the CRAC will be then computed as the solution of the optimization problem min mℓj,Tr,ℓj V(ℓ0,x,u,Tr) (7.12a) s.t. ˜xℓj+1 =h(xℓj, uℓj, Tr,ℓj,˜ Lkm,ℓj+1 ,˜ Wkm,ℓj+1 ) (7.12b) mℓj∈[1, M]∀j∈[0, Nc] (7.12c) uℓj= switch(mℓj)∀j∈[0, Nc] (7.12d) 15◦C≤Tr,ℓj≤25◦C∀j∈[0, Nc] (7.12e) Tℓj,(i)≤80◦C∀i∈M∀j∈[1, Np] (7.12f) mℓj≥mD,ℓj∀j∈[1, Np],(7.12g) where switch(·) represents the policy to select which servers are to be switched on. Here, a simple scheme is proposed for such policy. That is, for a given number of servers, the first mkwill be switched on. As it is customary in predictive controllers, the solution of (7.12) is applied in a receding horizon manner, meaning that only the first control actions and temperature set points are actually applied (i.e., uℓ0,(i)and Tr,ℓ0) while the remaining decision variables (uℓ1,(i), . . . , uℓNc−1,(i)and Tr,ℓ1, . . . , Tr,ℓNc−1) computed at time k are discarded. The optimization of (7.12) is then repeated at each sampling time so that the decision variables to be applied are computed using the real state of the data center at that sampling time. It should be noted that in this optimization problem some of the variables are integer (number of servers on) whereas other are real valued (the set-point temperatures). This, together with the complexity of the data center model motivates the use of specialized optimization algorithms to solve (7.12), such as a particle based optimization technique [180, 181] that will be exposed in the following section. 7.4 Particle based solvers for complex optimization problems In these techniques of iterative nature, a set of possible candidate solutions (called particles) are evaluated at each iteration and used to generate a new candidate solution set that may be closer to the solution of the optimization problem. Here, the evaluation of each particle will be done by means of simulations that will be used to assess the performance of each candidate solution. The main ingredients of the proposed technique are: •Particles are candidate solutions of the optimization problem. That is, a sequence of control actions over the control horizon. At higher problem 7.4. Particle based solvers for complex optimization problems 109 complexity (i.e., more servers and longer control horizons), more particles are needed in order to obtain a solution close to the optimal one [184]. •Weights. Particles have associated a weight that represents how good is a particle compared with the others. In this work, the weight is based on the performance cost of each particle computed by means of a computer simulation of the data center model. These simulations predict the evolution of the data center along the prediction horizon if the decision variables are those of the particles. Once the simulations for all particles are completed (this can be done in parallel), a scaling of the performance costs of the particles is made. Let Vmax and Vmin be the maximum and minimum costs attained on the set of particles under evaluation. If a given particle z≜ (m,Tr) has a performance cost V(z), then the weight of the particle zwill be σ(z)= 1 −V(z)−Vmin Vmax −Vmin ,(7.13) meaning that particles with higher costs will have lower weights and vice versa. This is very important for the resampling phase. •Feasibility checking. Those particles that do not satisfy the constraints in (7.12) will be assigned a zero weight so that they cannot be selected in the resampling stage. •Resampling is the step where a new generation of particles is created based on the performance (weights) at the previous iteration. This is done by means of the Kitagawa resampling algorithm with stratified scheme [185]. The aim of this method is to form groups of particles with good performance in different areas of the feasible set of solutions. In this way, the risk of getting stuck at a local minimum is reduced because the algorithm considers all particles and not only the best one. Particles with higher weights are more likely to be resampled at next iteration (i.e., selected to be included in the next particle set). •Perturbation. This process is inherently connected to the previous one. Once particles are resampled, it is necessary to ”move” them along the local search space in order to discover new feasible solutions with potentially lower costs. In this work, the perturbation is made by adding a Gaussian white noise. Figure 7.2 shows the main ideas of the resampling and perturbation steps in a 2 degree of freedom minimization example. At first, there exists a set of random-generated particles. As the iterations continue, particles move towards better possible solutions according to the costs obtained previously. In our case, the optimization variables are integer, thus this process ends rounding to the nearest integer. •Reseeding. The initial set of particles should be generated in a random manner. However, for the subsequent sample times, the reseeding can be 116 Chapter 7. Energy-efficient management of data centers 0 500 1000 1500 2000 2500 0 0.5 1 1.5 0 500 1000 1500 2000 2500 4 6 8 10 0 500 1000 1500 2000 2500 Time (s) 0 2 4 6 Figure 7.3: From top to bottom, Mean Queue Length, Mean Number of On Servers and Mean Number of Idle Servers controller. The top subplot shows the power consumption term in Watts (W). As it can be guessed, the controllers with larger number of active servers will have greater power consumption, due to the consumption of the servers themselves but also for the greater consumption of the CRAC unit (which following (7.9) is assumed proportional to that of the servers). The lower subplots are referred to the control efforts of the control actions. Here, it can be seen that the controllers with Nc= 5 tend to change more their control actions which leads to greater control efforts, but also facilitates a more tuned application of the control action that results in lower overall costs as seen in Figure 7.5. In this figure, the instantaneous total cost for each controller is shown. As expected, the more conservative controllers are the most expensive ones in terms of performance cost. 7.6.1 QoS violation rate Once all simulations are carried out, the bounds can be computed for each one of the four controllers. Table 7.2 shows the results for the QoS constraint. The final result is shown in the ρ+ ∆ρcolumns which are the upper bounds (with confidence 1 −δ) of the probability of not meeting the agreed service time Dfor the theorems 2.2 and 2.3 respectively. Lower numbers imply better probabilistic 7.6. Numerical results 117 0 500 1000 1500 2000 2500 1500 2000 2500 3000 3500 0 500 1000 1500 2000 2500 0 500 1000 0 500 1000 1500 2000 2500 Time (s) 0 200 400 Figure 7.4: From top to bottom: Mean Power Consumption, Mean Control Effort for booting servers and Mean Control Effort for changing the CRAC temperature reference guarantees, which as expected corresponds to the more conservative controllers (C1and C2). ρ ρ + ∆ρ(2.2) ρ+ ∆ρ(2.3) N C1S-Pbo, Nc= 1 0.0032 0.0050 0.0045 65960 C2S-Pbo, Nc= 5 0.0040 0.0060 0.0054 65926 C3Pbo, Nc= 1 0.0278 0.0328 0.0314 65708 C4Pbo, Nc= 5 0.0337 0.0392 0.0376 65678 Table 7.2: Violation rate of the proposed controllers for the QoS constraints with δ= 10−6 7.6.2 Thermal constraint violation rate Bounds have also been obtained for the temperature constraint, with table 7.3 summarizing the results. In this case, Nstands as the number of total time instants for all the simulations H(which were all of the same length), and the upper bound is on the probability of reaching a temperature higher than 80◦C. Again the more conservative a controller has the better probabilistic guarantees, but at the expense of a higher cost. 118 Chapter 7. Energy-efficient management of data centers 0 500 1000 1500 2000 2500 Time (s) 1800 2000 2200 2400 2600 2800 3000 3200 3400 3600 Figure 7.5: Instantaneous total cost ρ ρ + ∆ρ(2.2) ρ+ ∆ρ(2.3) N C1S-Pbo, Nc= 1 0.0000 0.0055 0.0058 2500 C2S-Pbo, Nc= 5 0.0000 0.0055 0.0058 2500 C3Pbo, Nc= 1 0.0001 0.0071 0.0067 2500 C4Pbo, Nc= 5 0.0000 0.0055 0.0058 2500 Table 7.3: Violation rate of the proposed controllers for the maximum temperature constraint with δ= 10−6 7.6.3 Parallel computation improvement In order to prove the speed-ups obtained with a parallel implementation, the algorithms have been implemented in the CPU and also in a CUDA capable GPU. The CPU version is coded completely in Matlab while the GPU version uses Matlab code for the serial operations and C code for the CUDA kernels. Table 7.4 shows the execution time of the controllers for an increasing number of particles. Cells with a hyphen mean that the time exceeded the sampling time. Finally, figure 7.6 shows the relative speed-up for different numbers of particles. S-Pbo GPU S-Pbo CPU Pbo GPU Pbo CPU 10 17.0931 23.4282 0.6974 1.0001 100 17.8046 - 0.7230 9.4192 1000 30.2881 - 1.2267 93.6546 10000 91.4658 - 3.7208 - Table 7.4: Mean computation time of the controllers in seconds. 7.6. Numerical results 119 0 100 200 300 400 500 600 700 800 900 1000 Number of Particles ( ) 0 10 20 30 40 50 60 70 80 Execution speed-up Figure 7.6: Speed-up obtained with the GPU computing for the controller execution. As it can be expected, at a higher number of particles, the execution time within the CPU grows strongly whereas the computation time within the GPU yields almost constant in comparison, which leads to the increasing slope shown in figure 7.6. It also should be noted that once reached a number of particles greater than 1000, the execution times of the CPU become unmanageable and the simulations are extremely costly. 7.6.4 Computation time analysis Also, an analysis of the computation time of the algorithm with respect to different number of servers has been done in order to study its behavior for the problem of large data centers. The results are shown in figure 7.7. A polynomial of order n= 2 is also added to the figures, confirming a quadratic complexity. 0 50 100 150 200 250 300 Number of Servers 0 50 100 150 200 250 Computation Time (s) Real Fitting 25 30 35 40 45 50 55 60 Number of Servers 20 40 60 80 100 120 140 160 180 200 Computation Time (s) Real Fitting Figure 7.7: Computation times of the algorithms for an increasing number of servers. Left: Algorithm without scenarios. Right: Algorithm with scenarios Thus, the algorithm is feasible from a computational point of view and it can be 120 Chapter 7. Energy-efficient management of data centers used with larger data centers. Nevertheless, for higher numbers of servers, besides adding more computational power, one could consider multiple controllers each one for each cold aisle or consider that the control actions handle clusters of servers instead of individual ones. 7.7 Conclusions This chapter presented a complex model of a data center and proposed a methodology to deal with the problem of minimizing the power consumption of the facilities while maintaining the quality of service at acceptable levels. The proposed MPC- inspired controller simulates the whole data center model in order to compute the optimal inputs to be applied to the system. For this reason, a particle-based algorithm was used to solve the optimization problem at each time instant k. The bounds proposed in chapter 6 were used to tune the hyperparameters of the controllers. 121 Chapter 8 Conclusions and future work This chapter briefly summarizes the contributions of this dissertation and presents some future lines of work in each field. 8.1 Contributions As it was stated in the introduction, the contributions of this thesis are related to three different fields: system identification, probabilistic forecasting and model predictive control. In what follows, the contributions of each chapter are stated: •Chapter 2 presented the notion of dissimilarity and the dissimilarity function to be used throughout this dissertation. It was shown that the proposed dissimilarity function has many advantageous properties. Also, by means of some forecasting examples, it was shown that the dissimilarity function can be used as a predictor, obtaining better results than those obtained with a neural network in the field of stock forecasting. •In chapter 3, the case of uncertainty predictions was tackled. Instead of providing just the expected value of the forecasting, a prediction region is obtained. First, a method oriented to obtain interval predictions of a univariate system or time series was proposed. This method was computationally expensive because it required the integration of the empirical conditioned probability density function for every time instant. In any case, the numerical examples showed that the performance was better in comparison to set-membership and quantile regression methods. Then, a method to bypass the numerical integration of the empirical conditioned probability density function was proposed. Besides this great advantage, which reduces considerably the execution time of the algorithm, multivariate prediction regions could be obtained. Also, it was shown that ellipsoidal approximations of such regions can be easily obtained. Finally, 122 Chapter 8. Conclusions and future work by means of a numerical example, it was shown that the proposed approach and the ellipsoidal approximations are smaller than those regions obtained with quantile regression, gaussian processes or inverse regression methods while fulfilling the probabilistic specifications. •Chapter 4 presented the state-space kriging method for autonomous systems. By using this technique, a model of the system can be obtained using only historical data of the process. Two different variations are proposed, one obtained by weighting the locality within the dissimilarity function and another one where the nonlinearity is modeled by means of kernel functions. Also, it was shown that the kalman filter can be used to improve the quality of the predictions. Finally, two numerical forecasting examples were presented, where it was shown that the proposed approaches obtain better results than other machine learning techniques. •On the other hand, chapter 5 presented the state-space kriging method for non-autonomous systems. Thanks to the modifications made to the constraints in the dissimilarity function, it is possible to include the effect of an external input. As the models obtained in this chapter are suitable for MPC, a MPC formulation using the K-SSK model was proposed. Finally, a simulation example and a real experiment were conducted in order to show the good performance of the controller. •In chapter 6, by means of the Chernoff’s bound [128, 159] and the results by Alamo et al. in [61], two different sharp bounds of the constraint violation rate of a certain controller were obtained. No assumptions on the probability distributions of the disturbances were made in order to prove the obtained results. It is only required to obtain i.i.d. samples, which is a standard assumption in the field. Finally, a simple MPC example was proposed to show how the proposed bounds should be computed. •Chapter 7 showcased a more complicated example where the proposed bounds can be used. The chosen system was a data center. This kind of facilities have an extremely large electrical consumption due to both IT equipment and infrastructure facilities. To tackle this problem, an energy-efficient management approach was proposed. Here, the bounds proposed in chapter 6 were used to tune the parameters of the controller. 8.2 Future work The future lines of work for each one of the parts of this dissertation can be found in the following: •As the obtained interval predictions and prediction regions of part I performed good enough in a forecasting example, it would be interesting to apply these uncertain predictions in a control scheme, i.e. developing new 8.2. Future work 123 stochastic MPC controllers whose prediction model is given by the predictors obtained in chapter 3. The most approachable way would be to synthesize a controller for a unidimensional system where the proposed interval predictor is used to obtain uncertain predictions of the system. These predictions would be used within the optimization problem to ensure the constraint satisfaction in a chanceconstrained manner. •The proposed MPC formulation in part II presented nominal stability and robust asymptotic stability, but lacked robust constraint satisfaction. For this reason, it would be interesting to improve the controller formulation to tackle this case. Also, the integration of a real time optimizer and the consideration of economic objectives in the controller is another pending task. On the other hand, the proposed SSK models hardly accept new data because this would mean that the state of the system is constantly increasing. It is necessary to find a good method that would allow the user to introduce new data to the SSK model without increasing the complexity. •When some of the proposed bounds in part III are used, it is possible to lose the recursive feasibility of the MPC controller unless the controller formulation is relaxed somehow. For that reason, a solution to this phenomenon should be proposed after a careful analysis of its origin. On the other hand, the i.i.d. assumption, although usual in the literature, is hard to be fulfilled because samples drawn from a dynamical system are inherently correlated. It would be interesting to develop methods which allow us to obtain samples that are as i.i.d. as possible. 124 Chapter 8. Conclusions and future work Bibliography [1] L. Ljung, System identification: Theory for the user. Pearson Education, 1998. [2] K. Ogata, Modern control engineering, vol. 5. Prentice hall Upper Saddle River, NJ, 2010. [3] P. Young, “Parameter estimation for continuous-time models—a survey,” Automatica, vol. 17, no. 1, pp. 23–39, 1981. [4] N. R. Kristensen, H. Madsen, and S. B. Jørgensen, “Parameter estimation in stochastic grey-box models,” Automatica, vol. 40, no. 2, pp. 225–237, 2004. [5] F. Ding, G. Liu, and X. P. Liu, “Parameter estimation with scarce measurements,” Automatica, vol. 47, no. 8, pp. 1646–1655, 2011. [6] P. Geladi and B. R. Kowalski, “Partial least-squares regression: a tutorial,” Analytica chimica acta, vol. 185, pp. 1–17, 1986. [7] P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” Journal of fluid mechanics, vol. 656, pp. 5–28, 2010. [8] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, Dynamic mode decomposition: data-driven modeling of complex systems. SIAM, 2016. [9] J. L. Proctor, S. L. Brunton, and J. N. Kutz, “Dynamic mode decomposition with control,” SIAM Journal on Applied Dynamical Systems, vol. 15, no. 1, pp. 142–161, 2016. [10] J. Schoukens and L. Ljung, “Nonlinear system identification: A useroriented road map,” IEEE Control Systems Magazine, vol. 39, no. 6, pp. 28– 99, 2019. [11] O. Nelles, “Nonlinear dynamic system identification,” in Nonlinear System Identification, pp. 547–577, Springer, 2001. [12] K. Liu, “Identification of linear time-varying systems,” Journal of Sound and Vibration, vol. 206, no. 4, pp. 487–505, 1997. 125 132 Bibliography [95] P. Krupa, I. Alvarado, D. Limon, and T. Alamo, “Implementation of model predictive control for tracking in embedded systems using a sparse extended ADMM algorithm,” IEEE Transactions on Control Systems Technology, 2021. [96] P. Krupa, D. Limon, and T. Alamo, “SPCIES: Suite of Predictive Controllers for Industrial Embedded Systems,” 2020. [97] P. O. Scokaert, D. Q. Mayne, and J. B. Rawlings, “Suboptimal model predictive control (feasibility implies stability),” IEEE Transactions on Automatic Control, vol. 44, no. 3, pp. 648–654, 1999. [98] G. Pannocchia, J. B. Rawlings, and S. J. Wright, “Conditions under which suboptimal nonlinear MPC is inherently robust,” Systems & Control Letters, vol. 60, no. 9, pp. 747–755, 2011. [99] M. M. Seron, J. A. De Dona, and G. C. Goodwin, “Global analytical model predictive control with input constraints,” in Proceedings of the 39th IEEE Conference on Decision and Control (Cat. No. 00CH37187), vol. 1, pp. 154– 159, IEEE, 2000. [100] A. Bemporad, M. Morari, V. Dua, and E. N. Pistikopoulos, “The explicit linear quadratic regulator for constrained systems,” Automatica, vol. 38, no. 1, pp. 3–20, 2002. [101] D. R. Ramirez and E. F. Camacho, “Piecewise affinity of min–max MPC with bounded additive uncertainties and a quadratic criterion,” Automatica, vol. 42, no. 2, pp. 295–302, 2006. [102] D. Mu˜noz de la Pe˜na, D. R. Ram´ırez, E. F. Camacho, and T. Alamo, “Explicit solution of min–max MPC with additive uncertainties and quadratic criterion,” Systems & control letters, vol. 55, no. 4, pp. 266–274, 2006. [103] E. Zafiriou, “On the closed-loop stability of constrained QDMC,” in 1991 American Control Conference, pp. 2367–2372, IEEE, 1991. [104] J. A. Primbs and V. Nevistic, “A framework for robustness analysis of constrained finite receding horizon control,” IEEE Transactions on Automatic Control, vol. 45, no. 10, pp. 1828–1838, 2000. [105] D. Limon, T. Alamo, D. M. Raimondo, D. Mu˜noz de la Pe˜na, J. M. Bravo, A. Ferramosca, and E. F. Camacho, “Input-to-state stability: a unifying framework for robust model predictive control,” in Nonlinear model predictive control, pp. 1–26, Springer, 2009. [106] G. Grimm, M. J. Messina, S. E. Tuna, and A. R. Teel, “Examples when nonlinear model predictive control is nonrobust,” Automatica, vol. 40, no. 10, pp. 1729–1738, 2004. Bibliography 133 [107] P. J. Campo and M. Morari, “Robust model predictive control,” in 1987 American control conference, pp. 1021–1026, IEEE, 1987. [108] D. R. Ramirez, T. Alamo, E. F. Camacho, and D. Mu˜noz de la Pe˜na, “Minmax MPC based on a computationally efficient upper bound of the worst case cost,” Journal of Process Control, vol. 16, no. 5, pp. 511–519, 2006. [109] P. O. Scokaert and D. Q. Mayne, “Min-max feedback model predictive control for constrained linear systems,” IEEE Transactions on Automatic control, vol. 43, no. 8, pp. 1136–1142, 1998. [110] D. Q. Mayne, E. C. Kerrigan, E. Van Wyk, and P. Falugi, “Tube-based robust nonlinear model predictive control,” International journal of robust and nonlinear control, vol. 21, no. 11, pp. 1341–1353, 2011. [111] D. Limon, I. Alvarado, T. Alamo, and E. F. Camacho, “Robust tube-based MPC for tracking of constrained linear systems with additive disturbances,” Journal of Process Control, vol. 20, no. 3, pp. 248–260, 2010. [112] R. Gonzalez, M. Fiacchini, T. Alamo, J. L. Guzman, and F. Rodr´ıguez, “Online robust tube-based MPC for time-varying systems: A practical approach,” International Journal of Control, vol. 84, no. 6, pp. 1157–1170, 2011. [113] D. Chatterjee, P. Hokayem, and J. Lygeros, “Stochastic receding horizon control with bounded control inputs: A vector space approach,” IEEE Transactions on Automatic Control, vol. 56, no. 11, pp. 2704–2710, 2011. [114] D. Chatterjee and J. Lygeros, “On stability and performance of stochastic predictive control techniques,” IEEE Transactions on Automatic Control, vol. 60, no. 2, pp. 509–514, 2014. [115] G. C. Calafiore and M. C. Campi, “The scenario approach to robust control design,” IEEE Transactions on automatic control, vol. 51, no. 5, pp. 742– 753, 2006. [116] M. C. Campi, S. Garatti, and M. Prandini, “Scenario optimization for MPC,” in Handbook of Model Predictive Control, pp. 445–463, Springer, 2019. [117] B. Karg, T. Alamo, and S. Lucia, “Probabilistic performance validation of deep learning-based robust nmpc controllers,” International Journal of Robust and Nonlinear Control, vol. 31, no. 18, pp. 8855–8876, 2021. [118] M. Alamir, “On probabilistic certification of combined cancer therapies using strongly uncertain models,” Journal of theoretical biology, vol. 384, pp. 59–69, 2015. 134 Bibliography [119] R. Amrit, J. B. Rawlings, and D. Angeli, “Economic optimization using model predictive control with a terminal cost,” Annual Reviews in Control, vol. 35, no. 2, pp. 178–186, 2011. [120] A. Ferramosca, J. B. Rawlings, D. Limon, and E. F. Camacho, “Economic MPC for a changing economic criterion,” in 49th IEEE Conference on Decision and Control (CDC), pp. 6131–6136, IEEE, 2010. [121] M. Diehl, R. Amrit, and J. B. Rawlings, “A Lyapunov function for economic optimizing model predictive control,” IEEE Transactions on Automatic Control, vol. 56, no. 3, pp. 703–707, 2010. [122] P. D. Christofides, R. Scattolini, D. Mu˜noz de la Pe˜na, and J. Liu, “Distributed model predictive control: A tutorial review and future research directions,” Computers & Chemical Engineering, vol. 51, pp. 21–41, 2013. [123] J. M. Maestre, R. R. Negenborn, et al.,Distributed model predictive control made easy, vol. 69. Springer, 2014. [124] D. Limon, I. Alvarado, T. Alamo, and E. F. Camacho, “MPC for tracking piecewise constant references for constrained linear systems,” Automatica, vol. 44, no. 9, pp. 2382–2387, 2008. [125] D. Limon, A. Ferramosca, I. Alvarado, and T. Alamo, “Nonlinear MPC for tracking piece-wise constant reference signals,” IEEE Transactions on Automatic Control, vol. 63, no. 11, pp. 3735–3750, 2018. [126] N. Cressie, “Kriging nonstationary data,” Journal of the American Statistical Association, vol. 81, no. 395, pp. 625–634, 1986. [127] J. P. Kleijnen, “Kriging metamodeling in simulation: A review,” European journal of operational research, vol. 192, no. 3, pp. 707–716, 2009. [128] H. Chernoff, “A measure of asymptotic efficiency for tests of a hypothesis based on the sum of observations,” The Annals of Mathematical Statistics, vol. 23, no. 4, pp. 493–507, 1952. [129] G. Alfonso, A. D. Carnerero, D. R. Ramirez, and T. Alamo, “Stock forecasting using local data,” IEEE Access, vol. 9, pp. 9334–9344, 2020. [130] A. D. Carnerero, D. R. Ramirez, and T. Alamo, “Probabilistic interval predictor based on dissimilarity functions,” IEEE Transactions on Automatic Control, 10.1109/TAC.2021.3136137, 2021. [131] A. D. Carnerero, D. R. Ramirez, and T. Alamo, “State-space kriging: A data-driven method to forecast nonlinear dynamical systems,” IEEE Control Systems Letters, vol. 6, pp. 2258 – 2263, 2022. [132] A. D. Carnerero, D. R. Ramirez, D. Limon, and T. Alamo, “Particle based optimization for predictive energy efficient data center management,” in Bibliography 135 2020 59th IEEE Conference on Decision and Control (CDC), pp. 2660– 2665, IEEE, 2020. [133] A. D. Carnerero, D. R. Ramirez, T. Alamo, and D. Limon, “Probabilistically certified management of data centers using predictive control,” IEEE Transactions on Automation Science and Engineering, 2021. [134] G. Alfonso, A. D. Carnerero, D. R. Ramirez, and T. Alamo, “Receding horizon optimization of large trade orders,” IEEE Access, vol. 9, pp. 63865– 63875, 2021. [135] A. Goshtasby, Image Registration. Principles, Tools and Methods. Springer, 2012. [136] S. T. Wierzcho´n and M. K lopotek, Modern algorithms of cluster analysis. Springer, 2018. [137] J. R. Salvador, D. R. Ramirez, T. Alamo, and D. Mu˜noz de la Pe˜na, “Offset free data driven control: application to a process control trainer,” IET Control Theory & Applications, vol. 13, no. 18, pp. 3096–3106, 2019. [138] J. R. Salvador, D. Mu˜noz de la Pe˜na, T. Alamo, and A. Bemporad, “Data-based predictive control via direct weight optimization,” IFAC- PapersOnLine, vol. 51, no. 20, pp. 356–361, 2018. [139] S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge University Press, 2004. [140] A. Beck, First-order methods in optimization. SIAM, 2017. [141] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sciences, vol. 2, no. 1, pp. 183–202, 2009. [142] Y. Nesterov, Lectures on convex optimization, vol. 137. Springer, 2018. [143] L. Ljung, System Identification. Prentice Hall, Upper Saddle River, NJ, 1999. [144] V. P. Upadhyay, S. Panwar, R. Merugu, and R. Panchariya, “Forecasting stock market movements using various kernel functions in support vector machine,” in Proceedings of the International Conference on Advances in Information Communication Technology —& Computing, AICTC ’16, (New York, NY, USA), Association for Computing Machinery, 2016. [145] M. Thomason, “The practitioner methods and tool,” Journal of Computational Intelligence in Finance, vol. 7, no. 3, pp. 36–45, 1999. [146] M. Mammarella, V. Mirasierra, M. Lorenzen, T. Alamo, and F. Dabbene, “Chance-constrained sets approximation: A probabilistic scaling approach,” Automatica, vol. 137, p. 110108, 2022. 136 Bibliography [147] T. Alamo, A. Cepeda, and D. Limon, “Improved computation of ellipsoidal invariant sets for saturated control systems,” in Proceedings of the 44th IEEE Conference on Decision and Control, pp. 6216–6221, IEEE, 2005. [148] K. P. Murphy, Machine learning: a probabilistic perspective. MIT press, 2012. [149] M. Erhard and H. Strauch, “Control of towing kites for seagoing vessels,” IEEE Transactions on Control Systems Technology, vol. 21, no. 5, pp. 1629– 1640, 2012. [150] N. Cristianini and J. Shawe-Taylor, An introduction to support vector machines and other kernel-based learning methods. Cambridge university press, 2000. [151] R. E. Kalman, “A New Approach to Linear Filtering and Prediction Problems,” Journal of Basic Engineering, vol. 82, no. 1, pp. 35–45, 1960. [152] H. K. Khalil, Nonlinear control, vol. 406. Pearson New York, 2015. [153] D. E. Seborg, T. F. Edgar, D. A. Mellichamp, and F. J. Doyle, Process dynamics and control. John Wiley & Sons, 2016. [154] G. Pannocchia and J. B. Rawlings, “Disturbance models for offset-free model-predictive control,” AIChE journal, vol. 49, no. 2, pp. 426–437, 2003. [155] U. Maeder, F. Borrelli, and M. Morari, “Linear offset-free Model Predictive Control,” Automatica, vol. 45, no. 10, pp. 2214–2222, 2009. [156] G. Pannocchia, “Offset-free tracking MPC: A tutorial review and comparison of different formulations,” in 2015 European Control Conference (ECC), pp. 527–532, 2015. [157] P. M. Oliveira and J. D. Hedengren, “An APMonitor temperature lab PID control experiment for undergraduate students,” in 2019 24th IEEE International Conference on Emerging Technologies and Factory Automation (ETFA), pp. 790–797, IEEE, 2019. [158] T. Alamo, R. Tempo, and E. F. Camacho, “Randomized strategies for probabilistic solutions of uncertain feasibility and optimization problems,” IEEE Transactions on Automatic Control, vol. 54, no. 11, pp. 2545–2559, 2009. [159] M. Okamoto, “Some inequalities relating to the partial sum of binomial probabilities,” Annals of the institute of Statistical Mathematics, vol. 10, no. 1, pp. 29–35, 1959. [160] H. Cram´er, “Les sommes et les fonctions de variables al´eatoires,” Actualit´es Scientifiques et Induestrielles. Conf´erences Internationales de Sciences. Paris Hermann, vol. 3, 1938. Bibliography 137 [161] A. Klenke, Probability theory: a comprehensive course. Springer Science & Business Media, 2013. [162] L. Parolini, B. Sinopoli, B. H. Krogh, and Z. Wang, “A cyber–physical systems approach to data center modeling and control for energy efficiency,” Proceedings of the IEEE, vol. 100, no. 1, pp. 254–268, 2011. [163] J. Scaramella, “Worldwide server power and cooling expense 2006-2010 forecast,” Market analysis, IDC Inc, 2006. [164] A. S. Andrae and T. Edler, “On global electricity usage of communication technology: trends to 2030,” Challenges, vol. 6, no. 1, pp. 117–157, 2015. [165] P. X. Gao, A. R. Curtis, B. Wong, and S. Keshav, “It’s not easy being green,” in Proceedings of the ACM SIGCOMM 2012 conference on Applications, technologies, architectures, and protocols for computer communication, pp. 211–222, 2012. [166] M. Dayarathna, Y. Wen, and R. Fan, “Data center energy consumption modeling: A survey,” IEEE Communications Surveys & Tutorials, vol. 18, no. 1, pp. 732–794, 2015. [167] W. Van Heddeghem, S. Lambert, B. Lannoo, D. Colle, M. Pickavet, and P. Demeester, “Trends in worldwide ICT electricity consumption from 2007 to 2012,” Computer Communications, vol. 50, pp. 64–76, 2014. [168] N. Lazic, C. Boutilier, T. Lu, E. Wong, B. Roy, M. Ryu, and G. Imwalle, “Data center cooling using model-predictive control,” in Advances in Neural Information Processing Systems, pp. 3814–3823, 2018. [169] H. Endo, S. Suzuki, H. Kodama, T. Hatanaka, H. Fukuda, and M. Fujita, “Development of predictive control system using just-in-time modeling and enthalpy-aware control in air conditioners for large-scale data center,” in 2018 18th International Conference on Control, Automation and Systems (ICCAS), pp. 1278–1283, IEEE, 2018. [170] M. Ogawa, H. Fukuda, H. Kodama, H. Endo, T. Sugimoto, T. Kasajima, and M. Kondo, “Development of a cooling control system for data centers utilizing indirect fresh air based on model predictive control,” in 2015 7th International Congress on Ultra Modern Telecommunications and Control Systems and Workshops (ICUMT), pp. 132–137, IEEE, 2015. [171] L. Fu, J. Wan, J. Yang, D. Cao, and G. Zhang, “Dynamic thermal and it resource management strategies for data center energy minimization,” Journal of Cloud Computing, vol. 6, no. 1, p. 25, 2017. [172] C. Bash and G. Forman, “Cool job allocation: Measuring the power savings of placing jobs at cooling-efficient locations in the data center.,” in USENIX Annual Technical Conference, vol. 138, pp. 140–149, 2007. 138 Bibliography [173] J. D. Moore, J. S. Chase, P. Ranganathan, and R. K. Sharma, “Making scheduling cool: Temperature-aware workload placement in data centers.,” in USENIX annual technical conference, pp. 61–75, 2005. [174] S. Li, H. Le, N. Pham, J. Heo, and T. Abdelzaher, “Joint optimization of computing and cooling energy: Analytic model and a machine room case study,” in 2012 IEEE 32nd International Conference on Distributed Computing Systems, pp. 396–405, 2012. [175] L. Parolini, B. Sinopoli, and B. H. Krogh, “Model predictive control of data centers in the smart grid scenario,” IFAC Proceedings Volumes, vol. 44, no. 1, pp. 10505–10510, 2011. [176] M. Ogura, J. Wan, and S. Kasahara, “Model predictive control for energyefficient operations of data centers with cold aisle containments,” IFAC- PapersOnLine, vol. 51, no. 20, pp. 209–214, 2018. [177] Q. Fang, J. Wang, and Q. Gong, “QoS-driven power management of data centers via model predictive control,” IEEE Transactions on Automation Science and Engineering, vol. 13, no. 4, pp. 1557–1566, 2016. [178] J. Yao, Z. Pan, and H. Zhang, “A distributed render farm system for animation production,” in Entertainment Computing – ICEC 2009 (S. Natkin and J. Dupire, eds.), pp. 264–269, Springer Berlin Heidelberg, 2009. [179] B. Bixby, “The Gurobi optimizer,” Transp. Research Part B, vol. 41, no. 2, pp. 159–178, 2007. [180] A. L. Visintini, W. Glover, J. Lygeros, and J. Maciejowski, “Monte Carlo optimization for conflict resolution in air traffic control,” IEEE Transactions on Intelligent Transportation Systems, vol. 7, no. 4, pp. 470–482, 2006. [181] N. Kantas, J. Maciejowski, and A. Lecchini-Visintini, “Sequential Monte Carlo for model predictive control,” in Nonlinear model predictive control, pp. 263–273, Springer, 2009. [182] G. Andreadis, L. Versluis, F. Mastenbroek, and A. Iosup, “A reference architecture for datacenter scheduling: design, validation, and experiments,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage, and Analysis, p. 37, 2018. [183] E. Chromy, T. Misuth, and M. Kavacky, “Erlang c formula and its use in the call centers,” Advances in Electrical and Electronic Engineering, vol. 9, no. 1, pp. 7–13, 2011. [184] J. P. De Villiers, S. Godsill, and S. Singh, “Particle predictive control,” Journal of Statistical Planning and Inference, vol. 141, no. 5, pp. 1753– 1763, 2011. Bibliography 139 [185] G. Kitagawa, “Monte Carlo filter and smoother for non-Gaussian nonlinear state space models,” Journal of computational and graphical statistics, vol. 5, no. 1, pp. 1–25, 1996. [186] J. Sanders and E. Kandrot, CUDA by example: an introduction to generalpurpose GPU programming, portable documents. Addison-Wesley Professional, 2010. [187] K. Karimi, N. G. Dickson, and F. Hamze, “A performance comparison of CUDA and OpenCL,” arXiv preprint arXiv:1005.2581, 2010. [188] S. Che, M. Boyer, J. Meng, D. Tarjan, J. W. Sheaffer, and K. Skadron, “A performance study of general-purpose applications on graphics processors using CUDA,” Journal of parallel and distributed computing, vol. 68, no. 10, pp. 1370–1380, 2008. [189] P. Bialas and A. Strzelecki, “Benchmarking the cost of thread divergence in CUDA,” in International Conference on Parallel Processing and Applied Mathematics, pp. 570–579, Springer, 2015.