程序代写代做 finance algorithm flex kernel C Bayesian rspa.royalsocietypublishing.org

程序代写 CS代考 / Algorithm算法代写代考, C语言代写, finance, kernel

rspa.royalsocietypublishing.org
Research
Cite this article: Kaiser E, Kutz JN, Brunton SL. 2018 Sparse identification of nonlinear dynamics for model predictive control in the low-data limit. Proc. R. Soc. A 474: 20180335. http://dx.doi.org/10.1098/rspa.2018.0335
Received: 23 May 2018 Accepted: 11 October 2018
Subject Areas:
applied mathematics, differential equations, mathematical modelling
Keywords:
model predictive control, nonlinear dynamics, sparse identification of nonlinear dynamics, system identification, control theory, machine learning
Author for correspondence:
E. Kaiser
e-mail: eurika@uw.edu
Sparse identification of
nonlinear dynamics for model
predictive control in the
low-data limit
E. Kaiser1, J. N. Kutz2 and S. L. Brunton1
1Department of Mechanical Engineering, and 2Department of Applied Mathematics, University of Washington, Seattle, WA, 98195
EK, 0000-0001-6049-0812
Data-driven discovery of dynamics via machine learning is pushing the frontiers of modelling and control efforts, providing a tremendous opportunity to extend the reach of model predictive control (MPC). However, many leading methods in machine learning, such as neural networks (NN), require large volumes of training data, may not be interpretable, do not easily include known constraints and symmetries, and may not generalize beyond the attractor where models are trained. These factors limit their use for the online identification of a model in the low-data limit, for example following an abrupt change to the system dynamics. In this work, we extend the recent sparse identification of nonlinear dynamics (SINDY) modelling procedure to include the effects of actuation and demonstrate the ability of these models to enhance the performance of MPC, based on limited, noisy data. SINDY models are parsimonious, identifying the fewest terms in the model needed to explain the data, making them interpretable and generalizable. We show that the resulting SINDY-MPC framework has higher performance, requires significantly less data, and is more computationally efficient and robust to noise than NN models, making it viable for online training and execution in response to rapid system changes. SINDY-MPC also shows improved performance over linear data-driven models, although linear models may provide a stopgap until enough data is available for SINDY. SINDY-MPC is demonstrated on a variety of dynamical systems with different challenges,
2018 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.

including the chaotic Lorenz system, a simple model for flight control of an F8 aircraft, and an 2 HIV model incorporating drug treatment.
1. Introduction
Data-fuelled modelling and control of complex systems is currently undergoing a revolution, driven by the confluence of big data, advanced algorithms in machine learning and modern computational hardware. Model-based control strategies, such as model predictive control (MPC), are ubiquitous, relying on accurate and efficient models that capture the relevant dynamics for a given objective. Increasingly, first principles models are giving way to data- driven approaches, for example in turbulence, epidemiology, neuroscience and finance [1]. Although these methods offer tremendous promise, there has been slow progress in distilling physical models of dynamic processes from data. Despite their undeniable success, many modern techniques in machine learning (e.g. neural networks (NN)) rely on access to massive datasets, have limited ability to generalize beyond the attractor where data is collected, and do not readily incorporate known physical constraints. The current challenges associated with data- driven discovery limit its use for real-time control of strongly nonlinear, high-dimensional, multi-scale systems, and prevent online recovery in response to abrupt changes in the dynamics. Fortunately, a new paradigm of sparse and parsimonious modelling is enabling interpretable models in the low-data limit. In this work, we extend the recent sparse identification of nonlinear dynamics (SINDy) framework [2] to identify models with actuation, and combine it with MPC for effective and interpretable data-driven, model-based control. We apply the proposed SINDY-MPC architecture to control several nonlinear systems and demonstrate improved control performance in the low-data limit, compared with other leading data-driven methods, including linear response models and NNs.
Model-based control techniques, such as MPC [3,4] and optimal control [5,6], are cornerstones of advanced process control, and are well-positioned to take advantage of the data-driven revolution. MPC is particularly ubiquitous in industrial applications, as it enables the control of strongly nonlinear systems with constraints, which are difficult to handle using traditional linear control approaches [7–11]. MPC benefits from simple and intuitive tuning and the ability to control a range of simple and complex phenomena, including systems with time delays, non-minimum phase dynamics, and instability. In addition, it is straightforward to incorporate known constraints and multiple operating conditions, it exhibits an intrinsic compensation for dead time, and it provides the flexibility to formulate and tailor a control objective. The major drawback of model-based control, such as MPC, lies in the development of a suitable model via existing system identification or model reduction techniques [12], which may require expensive and time-consuming data collection and computations.
Nearly all industrial applications of MPC rely on empirical models, and increasing plant complexity and tighter performance specifications require models with higher accuracy. There are many techniques to obtain data-driven models, including state-space models from the eigensystem realization algorithm (ERA) [13] and other subspace identification methods, Volterra series [14–16], autoregressive models [17] (e.g. ARX, ARMA, NARX and NARMAX [18] models), and NN models [19–22], to name only a few. These procedures all tend to yield black-box models, with limited interpretability, physical insights and ability to generalize. More recently, linear representations of nonlinear systems using extended dynamic mode decomposition [23] have been successfully paired with MPC [24,25]. Nonlinear models based on machine learning, such as NNs, are increasingly used due to advances in computing power, and recently deep reinforcement learning has been combined with MPC [26,27], yielding impressive results in the large-data limit. However, large volumes of data are often a luxury, and many systems must be identified and controlled with limited data, for example in response to abrupt changes. Current efforts are focused on rapid learning based on minimal data.
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

population dynamics (§3)
pitch … anglex2 rate x3
angle of attack
helper- dependent h
therapeutic strategies for infectious diseases (§6)
*
† l c1
d
† b2 c2q
helper- independent
x4
x1 1px5
surface
u
chaos control (§4)
automatic flight control (§5)
3
Figure 1. Applications of SINDY-MPC investigated in this work. (Online version in colour.)
When abrupt changes occur in the system, an effective controller must rapidly characterize and compensate for the new dynamics, leaving little time for discovery based on limited data. A second challenge is the ability of models to generalize beyond the training data, which is related to the ability to incorporate new information and quickly modify the model. Machine learning algorithms often suffer from overfitting and a lack of interpretability, although the application of these algorithms to physical systems offers a unique opportunity to incorporate known symmetries and constraints. These challenges point to the need for parsimonious and interpretable models [2,28,29] that may be characterized from limited data and in response to abrupt changes [30]. Whereas traditional methods require unrealistic amounts of training data, the recently proposed SINDY framework [2] relies on sparsity-promoting optimization to identify parsimonious models from limited data, resulting in interpretable models that avoid overfitting. It has also been shown recently [31] that it is possible to enforce known physics (e.g. constraints, conservation laws and symmetries) in the SINDY algorithm, improving stability and performance of models.
In this work, we combine SINDY with MPC for enhanced data-driven control of nonlinear systems in the low-data limit. First, we extend the SINDY architecture to identify interpretable models that include nonlinear dynamics and the effect of actuation. Next, we show the enhanced performance of SINDY-MPC compared with linear data-driven models and with NN models. The linear models are identified using dynamic mode decomposition with control (DMDc) [1,32], which is closely related to SINDY and traditional state-space modelling techniques such as ERA. SINDY-MPC is shown to have better prediction accuracy and control performance than NN models, especially for small and moderate amounts of noisy data. In addition, SINDY models are less expensive to train and execute than NN models, enabling real-time applications. SINDY- MPC also outperforms linear models for moderate amounts of data, although DMDc provides a working model in the extremely low-data limit for simple problems. Thus, in response to abrupt changes, a linear DMDc model may be used until a more accurate SINDY model is trained. We demonstrate the SINDY-MPC architecture on several systems of increasing complexity as illustrated in figure 1.
2. SINDY-MPC framework
The SINDY-MPC architecture combines the systematic data-driven discovery of dynamics with advanced model-based control to facilitate rapid model learning and control of strongly nonlinear
b
healthy cell
x CTL
c HIV
1 x3
1
b (1–hu)
precursor
infected cell
x2
† 1 a p2 †
†
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

SINDY-MPC model updating
e
trust region: ||yj – yˆj||2 0 and R􏰋u > 0 are positive definite and positive semi-definite, respectively. Note that the model prediction xˆk, which is forecast, may differ from the true measured state xk.
The dynamics are given by the identified SINDYc model, e.g. x ̇=F(x,u)=ΞΘT(x,u); Fˆ represents a discrete-time or discretized SINDYc model. While the model and the control law may be learned simultaneously, we adopt a two-stage process, where the model is first learned from data and then used in the control optimization with MPC. A joint optimization of the model and the control law may be challenging, as the particular control action depends on the model. However, it may be possible to develop a streaming algorithm to adapt the model to abrupt system changes [30], iterating between model identification and control optimization.
MPC is one of the most powerful model-based control techniques due to the flexibility in the formulation of the objective functional, the ability to add constraints, and extensions to nonlinear systems. The most challenging aspect of MPC involves the identification of a dynamical model that accurately and efficiently represents the system behaviour when control is applied. If the model is linear, minimization of a quadratic cost functional subject to linear constraints results in a tractable convex problem. Nonlinear models may yield significant improvements; however, they render MPC a nonlinear program, which can be expensive to solve, making it particularly challenging for real-time control. Conditions on the well-posedness of the problem and existence and uniqueness of the solution of the nonlinear optimization problem are, e.g. provided in [60]. Fortunately, improvements in computing power and advanced algorithms are increasingly enabling nonlinear MPC for real-time applications.
3. A simple model for population dynamics
We first demonstrate the SINDY-MPC architecture on the Lotka–Volterra system, a two- dimensional, weakly nonlinear dynamical system, describing the interaction between two competing populations. These dynamics may represent two species in biological systems, competition in stock markets [61], and can be modified to study the spread of infectious diseases [62]. We will consider more sophisticated examples in the following sections.
The dynamics of the prey and predator populations, x1 and x2, respectively, are given by
x ̇1 = ax1 − bx1x2 (3.1a)
and
x ̇2 = −cx2 + dx1x2 + u, (3.1b)
where the constant parameters a = 0.5, b = 0.025, c = 0.5 and d = 0.005 represent the growth/death
rates, the effect of predation on the prey population, and the growth of predators based on the size
of the prey population. The unforced system exhibits a limit cycle behaviour, where the predator
lags the prey, and a critical point xcrit = (g/d a/b)T, where the population sizes of both species are
in balance. The control objective is to stabilize this fixed point. Here, the timestep 􏰋t = 0.1 of the
system and the model are equal, the weight matrices are Q = ( 1 0 ) and Ru = R􏰋u = 0.5, and the 01
actuation input is limited to u ∈ [−20, 20]. The control and prediction horizons are mp = mc = 10 unless otherwise noted. We apply an additional constraint on u, so that x2 does not decrease below 10, to enforce a minimum population size required for recovery.
To assess the performance and capabilities of the SINDY-MPC architecture, SINDYc is compared with two representative data-driven models: dynamic mode decomposition with control (DMDc) and a multilayer NN, which can represent any continuous function under mild conditions [63]. The results are displayed in figure 5. The first 100 time units are used to train the models with a phase-shifted sum of sinusoids as input, a so-called Schroeder sweep [64],
9
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

(a)
training
prediction control DMDc
SINDYc
NN
10
100
0
(c)
(b)
6
4 2 0
200 250 time
execution time
0.39 s 0.31 s
19 min
DMDc SINDYc NN
300
0
100 predator
200 300 prey control
time
Figure 5. Prediction and control performance for the Lotka–Volterra system: (a) time series of states and input during training, validation and control stage, (b) cumulative cost and (c) execution time of the MPC optimization procedure. (Online version in colour.)
after which the predictive capabilities of these models are validated using sinusoidal forcing with u(t) = (2 sin(t) sin(t/10))2 on the next 100 time units. Different actuation inputs are used during the training and validation stages to assess the models’ ability to generalize. Thereafter, MPC is applied for 100 time units using a prediction and control horizon of mp = mc = 5. SINDYc shows the best prediction and control performance, followed by DMDc and the NN (due to its steady- state error). The NN has 1 hidden layer with 10 neurons, which is the best trade-off between model complexity and accuracy; increasing the number of neurons or layers has little impact on the prediction performance. Further, hyperbolic tangent sigmoid activation functions are employed. It is first trained as a feed-forward network using the Levenberg–Marquardt algorithm and then closed. If the data are corrupted by noise, a Bayesian regularization is employed, which requires more training time but improves robustness. While the NN exhibits a similar control performance, the execution time of SINDYc is 37 times faster, which is particularly critical in real- time applications. For a fair comparison, all methods are compared using the same optimization routine based on interior-point methods via Matlab’s fmincon. Thus, it would be possible to reduce the time for the linear system further.
In practice, measurements are generally affected by noise. We examine the robustness of these models for increasing noise corruption of the state measurements, i.e. y = x + n where n ∈ N(0,σ2) with standard deviation σ. Cross-validated prediction performance for different noise magnitudes η = σ/ max(std(xi)) ∈ (0.01, 0.5), where std denotes standard deviation, is displayed in figure 7a,b. As expected, the performance of all models decreases with increasing noise magnitude. SINDYc generally outperforms DMDc and NN models, exhibiting a slower decline in performance for low and moderate noise levels. Sparse regression is known to improve robustness to noise and prevent overfitting. The large fluctuation in the NN performance are due to its strong dependency on the initial network weights.
The amount of data required to train an accurate model is particularly crucial in real-time applications, where abrupt changes or actuation may render the model invalid and rapid model updates are necessary. Figure 6a–c shows the average relative prediction error on 100 time units used for validation, and the training time for increasing lengths of training data. The effect of the
population size
model complexity
cost (×104)
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

(a)
(d)
107 106 105 104
11
100 50 0
predator prey
DMDc SINDYc NN
h=0
102 103 training length m
(b)
1
10–1 10–2
(c)
1
10–1 10–2 10–3
0
10
10
(e)
80
20
10
9×103 8×103 7 × 103
200
100
50 100 time t/training length m × 10
205
JT=20 210 215 220
time
Figure 6. Cross-validated prediction error and control performance for increasing length of training data: (a) training time series, (b) average relative prediction error, (c) model training time in seconds, (d) terminal cumulative cost over 20 time units and (e) time series of states and cost of the best model for each model type (mDMDc = 20, mSINDYc = 85, mNN = 103 ). From m = 14 onwards, SINDYc yields highly performing models in MPC, outperforming all other models. Outside the shaded region, models perform significantly worse or even diverge. (Online version in colour.)
training length on the control performance (evaluated over 20-time units) is shown in figure 6d,e. For small amounts of data, the sparsity-promoting parameter λ in SINDYc is reduced by a factor of 10 until a non-zero entry appears. In the low-data limit, a highly predictive SINDYc model can be learned, discovering the true governing equations within machine precision. Significantly larger amounts of data are required to train an accurate NN model, although with enough data it outperforms DMDc. DMDc models may be useful in the extremely low-data limit, before enough data is available to characterize a SINDYc model. The training times of SINDYc and DMDc models increase slightly with the amount of data, but they require about two orders of magnitude less time than NN models. SINDYc’s intrinsic robustness to overfitting renders all models from mtrain = 14 on as having the best control performance compared with the overall best performing DMDc and NN models. By contrast, DMDc shows a slight decrease in performance due to overfitting and the NN’s dependency on the initial network weights detrimentally affects its performance. It is interesting to note that the control performance is generally less sensitive than the long-term prediction performance shown in figure 6b,c. Even a model with moderately low predictive accuracy may perform well in MPC.
In figure 7c,d, we show the same analysis but with noise-corrupted training data. We assume no noise corruption during the control stage. For each training length, the best model out of 50 noise realizations is tested for control. DMDc and SINDYc models both require slightly more data to achieve a similar performance as without noise. Note that NN models perform significantly worse when trained on noise-corrupted data.
4. Chaotic Lorenz system
In this section, we demonstrate the SINDY-MPC architecture on the chaotic Lorenz system, a prototypical example of chaos in dynamical systems. The Lorenz system represents the
training time (s) avg. rel. error population size
cost predator prey
cost JT = 20
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

(a)
0 (b) 1.0
0.5
0 0.01 0.1
h=0.05
102 103 length of training data
DMDc
SINDYc
NN
(c)
107 106 105 104
12
1000
500
10
0.2
0.3
0.4
0.5
(d)
100
80
20 10
9×103 8×103 7×103200
205 210 time
JT=20 215 220
h
Figure 7. Cross-validated prediction and control performance for increasing measurement noise: (a) mean squared and (b) average relative prediction error, (c) terminal cumulative cost over 20 time units and (d) time series of states and cost of the best models (mDMDc = 12, mSINDYc = 1250, mNN = 65). The control performance is shown for increasing length of noise-corrupted training data with η = 0.5. From m = 200 onwards, SINDYc yields highly performing models, outperforming all other models. Statistics are shown for 50 noise realizations each. Above the shaded region, models in most realizations do not have any predictive power. (Online version in colour.)
Rayleigh–Bénard convection in fluid dynamics as proposed by Lorenz [65], but has also been associated with lasers, dynamos and chemical reaction systems. The Lorenz dynamics are given by
x ̇1 =σ(x2 −x1)+u
x ̇2 = x1(ρ − x3) − x2 and x ̇3 = x1x2 − βx3
(4.1a) (4.1b) (4.1c)
with system parameters σ = 10, β = 8/3, ρ = 28, and control input u affecting only the first state. A typical trajectory oscillates alternately around the two weakly unstable fixed points (±√72,±√72,27)T. The chaotic motion of the system implies a strong sensitivity to initial conditions, i.e. small uncertainties in the state will grow exponentially with time. This represents a particularly challenging problem for model identification and subsequent control, as measurement and model uncertainty both lead to long-time forecast error.
The control objective is to stabilize one of these fixed points. In general, the timestep of the model is chosen to balance the control horizon, the length of the sequence of control inputs to be optimized and prediction accuracy. Here, the system timestep is 􏰋tsys = 0.001 and the model timestep is 􏰋tmodel = 0.01. The control input is determined every 10 system timesteps and then held constant. The weight matrices are Q = I3, where In denotes a n × n identity matrix, Ru = R􏰋u = 0.001, and the actuation input is limited to u ∈ [−50, 50]. The control and prediction horizon is mp = mc = 10 and the sparsity-promoting parameter in SINDYc is λ = 0.1, unless otherwise noted. For all cases, we assume access to full-state information.
average rel. error mean squared error
cost predator prey cost JT = 20
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

(a)
training
prediction control DMDc
SINDYc
NN
13
(b)
3
2 1 0
DMDc SINDYc NN
40 20 0 –20 –40
(c)
20 21 22 23 24 25 time
execution time (in min)
1.9 1.6
35
0 5 10 15 20 25 time
x1 x2 x3
control
Figure 8. Prediction and control performance for the chaotic Lorenz system: (a) time series of the states and input (shifted to −25 and scaled by 10 to improve readability) during training, validation and control stage, (b) cumulative cost and (c) execution time of the MPC optimization. (Online version in colour.)
We compare the prediction and control performance of the SINDYc model with DMDc and NN models. DMDc is trained to model the deviation from the goal state by constructing the regression model based on data from which the goal state has been subtracted. A less naive approach would partition the trajectory into two bins, e.g. based on negative and positive values of x1, and estimate two models for each goal state separately. The NN consists of 1 hidden layer with 10 neurons and employs hyperbolic tangent sigmoid activation functions. Cross-validated prediction and control performance for the Lorenz system are displayed in figure 8. The first 10 time units are used to train with a Schroeder sweep, after which the models are validated on the next 10 time units using a sinusoidally based high-frequency forcing, u(t) = (5 sin(30t))3. MPC is then applied for the last 5 time units. SINDYc exhibits the best prediction and control performance. The NN exhibits comparable control performance, although the prediction horizon is considerably shorter. Surprisingly, DMDc is able to stabilize the fixed point, despite poor predictions based on a linear model. As the predictive capability of DMDc is poor, we will not present DMDc results in the following, but instead compare SINDYc and the NN. As in the previous example, while the NN exhibits similar control performance, the control execution of SINDYc is 21 times faster.
Figure 9 examines the cross-validated prediction performance of SINDYc and NN models based on measurements with increasing noise magnitude η = σ/ max(std(xi)) ∈ {0.01, 0.1, 0.25}. The performance of both models decreases with increasing noise level, although SINDYc generally outperforms the NN. Unlike the Lotka–Volterra model, the average relative error is misleading in this case. With increasing noise magnitude the NN converges to a fixed point, having no predictive power, while SINDYc still exhibits the correct statistics beyond the prediction horizon; however, a phase drift leads to a larger average relative error. This is shown in figure 9b with the median (thick line) and the 25–75 percentile region (shaded area) of the prediction for three different noise levels. Thus, a better metric for prediction performance is the prediction horizon itself (figure 9a). The prediction horizon is estimated as the time instant when the error ball is larger than a radius of ε = 3, i.e. a model is considered predictive if
xi
model complexity
cost (×104)
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

(a)
10
10–1 10–3
NN
SINDYc
0.01 0.1 0.2 0.3 0.4 0.5
h
x3 x2 x1
(b)
h = 0.01
h = 0.1
h = 0.25
14
30
time
40
Figure 9. Cross-validated prediction performance for increasing measurement noise: (a) prediction horizon in time units and (b) time series with 50 (median as thick coloured line) and 25–75 (coloured shaded region) percentiles. Statistics are shown for 50 noise realizations each and different noise levels η ∈ {0.01, 0.1, 0.25}. (Online version in colour.)
􏰌􏰆3i=1(xi − xˆi)2 < ε. This corresponds to roughly 10% error per state variable, considering that the order of magnitude of each state is approximately O(101); this error radius correlates well with the visible divergence of the true and predicted state in figure 9b. For low and moderate noise levels, SINDYc robustly predicts the state with high accuracy. Even for η = 0.25, the 1-period prediction is sufficiently long for a successful stabilization with MPC as we consider a comparably short prediction horizon of Tp = 0.1. The effect of the amount of training data on the prediction and control performance is examined in figure 10, respectively. In figure 10a–d, we show the average relative error evaluated on the prediction over the next 10 time units, the prediction horizon, and the required training time in seconds for increasing length of noise-free training data. For a relatively small amount of data, SINDYc rapidly outperforms the NN model with a prediction horizon of 2.5 time units and a significantly smaller error. For a sufficiently large amount of data, SINDYc and the NN result in comparable predictions. However, SINDYc yields highly predictive models that can be rapidly trained in low and moderate data regimes. Models trained on weakly noise-corrupted measurements, η = 0.05, are tested in MPC. For each length of training data, 50 noise realizations are performed and the most predictive model is selected for evaluation in MPC (figure 10e,f). Outside the shaded regions, models are generally not predictive or might even diverge. In the noise-corrupted case, it is clear that SINDYc models generally have better control performance than NN models. For a sufficiently large amount of training data, NNs can have comparable performance to SINDYc models, although they show a sensitive dependence on the initial choice of the network weights. The control results of the NN are significantly better here than for the Lotka–Volterra model due to the intrinsic system properties. In chaotic systems, a long enough trajectory will come arbitrarily close to every point on the attractor; thus, measurements of the Lorenz system are in some sense richer than those of the Lotka–Volterra system. A surprising result is that a nearly optimal SINDYc model can be trained on just eight noisy measurements (compare figure 10e,f ). 5. Tracking for the F-8 crusader In this section, we consider an automatic flight control system of the F-8 aircraft at an altitude of 30 000 ft (9000 m) and Mach = 0.85 [66–68]. The control objective is to track a specific trajectory of increasing noise prediction horizon NN SINDYc rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ................................................... 60 40 20 0 –20 (b) (c) (d) (a) x1 (e) x2 x3 NN 15 106 104 102 1 10 1 10–1 10–2 10 10–1 10–3 SINDYc 105 (f) 0 –5 –10 0 –10 –20 40 20 0 3 2 1 h = 0.05 102 103 length of training data 10 0 5 10 time t / training lenght m × 100 JT=3 20 21 22 23 time Figure 10. Cross-validated prediction and control performance for increasing length of training data (without noise): (a) time series of the training data, (b) average relative prediction error, (c) prediction horizon, (d) training time in seconds, (e) terminal cumulative cost over 3 time units and (f ) time series of states and cost of the best model for each model type (mSINDYc = 38, mNN = 40). Note that from m = 400 onwards, SINDY identifies the best performing models. (Online version in colour.) the angle of attack. The aircraft dynamics [66] are given by x ̇ 1 = −0.877x1 + x3 − 0.088x1 x3 + 0.47x21 − 0.019x2 − x21 x3 + 3.846x31 − 0.215u + 0.28x21u + 0.47x1u2 + 0.63u3 x ̇2 =x3 and x ̇ 3 = −4.208x1 − 0.396x3 − 0.47x21 − 3.564x31 − 20.967u + 6.265x21 u + 46x1 u2 + 61.1u3 (5.1a) (5.1b) (5.1c) where x1 is the angle of attack (rad), x2 is the pitch angle (rad), x3 is the pitch rate (rad s−1 ) and u is the control input representing the tail deflection angle (rad). The system is non-affine in the states and the control input rendering it strongly nonlinear. The commanded angle of attack to be tracked [68] is given by 􏰃 0.5 1 􏰄 r(t)=0.4 −1+eˆt−0.8 + 1+eˆt−3 −0.4 (5.2) with ˆt = t/0.1. We assume that the output, over which the performance is optimized, is y = x1. The timestep of the system is 􏰋t = 0.001 and the timestep of the model is 􏰋tM = 0.01. The control input is determined using SINDY-MPC every 10 system timesteps over which the applied control is then kept constant. The weight matrices are Q = 25, Ru = R􏰋u = 0.05, the actuation input rate is limited to 􏰋u∈[−0.3,0.5], and the constraint for the angle of attack is y∈[−0.2,0.4]. The control and prediction horizon is mp = mc = 13 and the sparsity-promoting parameter in SINDYc rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ................................................... training time (s) prediction avg. rel. error states xi horizon cost (×104) x3 x2 x1 cost JT = 3 (a) training DMDc 1 0 –1 SINDYc –2 4 2 0 –2 –4 prediction (b) (c) 0.4 0.2 0 –0.2 –0.4 (d) 1.5 0.5 0 150 152 16 NN 0 50 100 150 time x1 x2 x3 control 0.2 0.1 0 –0.1 –0.2 DMDc SINDYc NN r Figure 11. Prediction and control performance for F8 aircraft: (a) time series of states and input during training, validation and control, (b) angle of attack y = x1 with reference r, (c) control input and (d) cumulative cost during control. (Online version in colour.) is λ = (10−4, 10−2, 10−2), where λi is used to identify the terms for xi. The NN has two hidden layers each with 15 neurons. Access to full-state information is assumed for these models. Results assessing prediction and control performance of SINDYc compared with DMDc and a NN model are displayed in figure 11. Similar to the Lotka–Volterra system, the NN requires more and richer training data, i.e. a better exploration of the system behaviour, to perform sufficiently well. Thus, 250 short trajectories each consisting of 1000 snapshots (25 × 104 instances in total) with varying input signals are used to train the NN; a subset of 20 trajectories is displayed in figure 11a(bottom). By contrast, SINDYc and DMDc perform similarly well if trained on much less data (104 instances of a single trajectory). Moreover, SINDYc learns from few measurements the true relationship between the variables, even though only limited system behaviour has been observed, resulting in increased performance. 6. Optimal therapy for pathogenic attacks Optimizing drug therapy is critical for inhibiting diseases such as cancer and viral infections. Here, we consider treatment of infections with the human immunodeficiency virus (HIV), a pathogen that infects T-helper CD4+ cells of the immune system and can cause acquired immune deficiency syndrome (AIDS). Identifying the underlying infection mechanism, the response of the immune system, and the interactions with drugs targeting different components in this system is critical for developing and optimizing therapeutic strategies. Various models have been proposed to study the interaction between HIV and CD4+ cells; we refer to a recent review [69]. Optimal treatment aims to decrease virus mutations, complications from administered drugs, medical costs, and to strengthen the immune system. We consider a system [70] that incorporates time 154 156 cost model complexity control input angle of attack xi xi rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ................................................... infections with HIV, the cytotoxic lymphocyte (CTL) response of the immune system, and therapeutic interventions via a highly active anti-retroviral therapy (HAART), i.e. a combination of drugs that affect the replication rate of HIV and support the immune system. This is based on a more general and complex system, that can be simplified under certain conditions [71]: 17 x ̇1 = λ − dx1 − β(1 − ηu)x1x2 x ̇2 = β(1 − ηu)x1x2 − ax2 − p1x4x2 − p2x5x2 x ̇3 = c2x1x2x3 − c2qx2x3 − b2x3 x ̇4 = c1x2x4 − b1x4 and x ̇5 = c2qx2x3 − hx5 (6.1a) (6.1b) (6.1c) (6.1d) (6.1e) with parameters λ=1, d=0.1, β =1, a=0.2, p1 =1, p2 =1, c1 =0.03, c2 =0.06, b1 =0.1, b2 = 0.01, q=0.5, h=0.1 and η=0.9799 (units typically in mm−3d−1). Here, the states describe concentrations of healthy CD4+ T-cells, x1, HIV-infected CD4+ T-cells, x2, CTL precursors (memory CTL), x3, helper-independent CTL, x4 and helper-dependent CTL, x5. For a detailed discussion of the system (6.1) we refer to [70,71]. The parameter η represents the effectiveness of the HAART therapy applied via u. For the considered parameters and in the absence of control (u ≡ 0), the system exhibits two stable fixed points: a progressive infection leading to AIDS, xA , and the recovery from a successful immune response, xB. The later steady state is given by xB = λ , xB = c2(λ−dq)−b2β − 1 d+βxB 2 2 xB = hxB5 , 3 c2qxB [c2(λ−dq)−b2β]2 −4βc2qdb2 2βc2q (6.2a) (6.2b) 􏰌 and xB =0, 4 xB = xB2c2(βq−a)+b2β, 5 c2p2xB 22 and exists if [c2(λ − dq) − b2β]2 − 4βc2qdb2 ≥ 0. The region of attraction (ROA) to this fixed point is limited and only established if the infectivity of the virus is small such that β < c1[c2b2(λ − qd) − b2c1d]/b1(c2b1q + b2c1), which can be achieved by applying a HAART therapy (u=1)withhighefficacy(η≈1).Thisstatemovesasafunctionofβeff =β(1−ηu)whenu>0and its ROA changes and does not necessarily overlap with the ROA in the absence of drug treatment (u = 0), i.e. dependent on the initial condition and the applied control the system will converge to a different steady state. A non-trivial control strategy is required that switches between treatment and no treatment to establish a successful immune response, and hence to approach xB. By contrast, when treatment is applied continuously for a sufficiently long amount of time such that the fixed points are approached and then terminated, the system will converge to a progressive infection, xA, even if a successful immune response had been established.
The cost functional to be optimized is given by
􏰑T 0
where xˆ 1 = xB1 and xˆ 3 = xB3 [70] taking into account the healthy cells, the immune system, and the cost of treatment. The control input u is bounded by 0 ≤ u ≤ 1 with efficacy of η = 0.9799. An additional constraint is added to the control that renders all cell concentrations non-negative, i.e. xi ≥0∀i. The time step is 􏰋tM =2h for the model and is 􏰋t=1/24day=1h for the simulated system. The control performance is evaluated over 50 weeks. The prediction and control horizon for the MPC optimization are both mp = mc = 24, i.e. over 2 days (from mp􏰋tM). We assume a more realistic situation, where the state is measured once a week, and the treatment is then kept constant over the following week. The training data consists of samples collected over 200 days (≈ 30 weeks) with a discrete control input, as was applied for the validation data in figure 12 (bottom). By contrast, the training data for the NN consists of ensemble data of 32 different
J =
(x1(t) − xˆ1) + (x3(t) − xˆ3) + |u(t)| dt, (6.3)
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

(a) x1
6
4 x2 2 0 –2
40 x3 20 0 1.0
x4 0.5
prediction
(b)
xi
xi
control SINDYc
DMDc
eDMDc
control 18 PI-SINDYc
delay-DMDc
NN
10 5 0
0.8 0.4 0
0.8
0
4 20 0
0.4
x5
–2
105 10–5 10–15
u 0.5
0 5 10
time (weeks)
0.8
xi
0.4 0
10 6 2
x1 x2 x3 xB1 xB2 xB3
DMDc delay-DMDc eDMDc SINDYc PI-SINDYc NN
Figure 12. Prediction and control performance for various models of the HIV model: (a) Predicted and true (solid black) time series of the forced dynamics and (b) control results (normalized for better legibility) with optimized actuation input (shaded background) and reference values. Only SINDYc and, with slightly less performance, PI-SINDYc and delay-DMDc are capable of driving the system to the desired steady state. (Online version in colour.)
trajectories. In both cases, the control is a random sequence of values that are kept constant over random durations of time 􏰋T ∈ [5 h, 10 days].
We consider a SINDYc model with (1) full-state information (SINDYc) and (2) partial information based on a subset of the variables (PI-SINDYc). The latter case demonstrates the situation when only a few states can be measured, which is generally more realistic. For the identification of the SINDYc models, it is important to normalize first the features in the library, as the coefficients of the active terms spread over several orders of magnitude. In both cases, a polynomial order of three is used for the library. The results are compared with various linear models: (1) DMDc on the full state (DMDc), (2) DMDc on delay coordinates of the full state (Delay- DMDc) and (3) DMDc on a set of nonlinear observables (extended DMD with control, eDMDc). In addition, a NN model on the full state (NN) is trained. An overview of these models and their parameters is provided in table 1. The identified parameters of the SINDYc and PI-SINDYc
0
10 20 30 40 50
time (weeks)
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………
cost (×109)
error

Table 1. Model parameters for the HIV system.
model state settings # unknowns
SINDYc y = x polynomial basis (order r = 3), 84 × 5 = 420 λ = (10, 3.1, 3, 0.1, 0.5)
………………………………………………………………………………………………………………………………………………………………………………….
PI-SINDYc y = (x1, x2, x3) polynomial basis (order r = 3), λ = (10, 30, 3) 35 × 3 = 105 ………………………………………………………………………………………………………………………………………………………………………………….
DMDc y = x − xB deviation from reference state 52 = 25 ………………………………………………………………………………………………………………………………………………………………………………….
Delay-DMDc y = (x(t) − xB , x(t − τ ) − deviation from reference state, 10 time delay (10 × 5)2 = 2500 xB,…, x(t − 9τ) − xB) coordinates of the full-state and of the control
input ………………………………………………………………………………………………………………………………………………………………………………….
eDMDc y = Θ(x, r) order r = 3 of polynomial basis (without 832 = 6889 constant term)
…………………………………………………………………………………………………………………………………………………………………………………. NN y = x 1 hidden layer with 5 neurons and linear 43
activation functions, data is log-transformed and mapped to [−1, 1] to compensate for skewness and different range
………………………………………………………………………………………………………………………………………………………………………………….
19
models are displayed below (6.4) are:
⎡0.9995 0 0 0
⎢-0.0999 0 0 0
⎢ 0 -0.1991 0 0
⎢ 0 0 -0.0100 0
⎢ 0 0 0 -0.1000
⎢ 0 0 0 0 -0.1000⎥ ⎢ 0 0 0 ⎥ ⎢ -0.9990 0.9981 0 0 0 ⎥ ⎢ -0.9990 0.9424 0 ⎥
⎢ 0
⎢ 0
⎢ 0
⎢ 0
⎢⎣ 0
0 -0.0299 0 0.0300 ⎥ ⎢
0 -0.0573 -0.0299 ⎥ 0 00⎥ 0 0 0 ⎥ 0 0.0069 0.0600 ⎥ 0 0 0 ⎥⎦ 0 0 0
0.9763 -0.7507 0
-0.9982 0 0.0300 -0.9990 0 0
0 0.0600 0 00 0 00 0
0 ⎥ ⎢ 0 ⎥ ⎢ 0 ⎥ ⎢ 0 ⎥⎦ ⎢⎣ 0
0
0.9763 -0.9757 0 0 0
􏰍 􏰎􏰏
Ξ SINDYc
􏰐􏰍 􏰎􏰏 􏰐
Ξ PI−SINDYc
0⎤⎡0.999500⎤ 0 ⎥ ⎢-0.0999 0 0 ⎥ 0 ⎥⎢ 0 -0.8985 0 ⎥ 0 ⎥ ⎢ 0 0 -0.0100⎥ 0 ⎥ ⎢ 0 0 0 ⎥
1 x1 x2 x3 x4 x5 x1x2 x2x3 x2 x4 x2 x5 x1x2x3 x1x2x4 x1x2x5 x1x2u
(6.4)
Only the non-zero parameters are shown, and their error is O(10−3)−O(10−6) for SINDYc. The error in the parameters decreases with increased time resolution. Here, a coarse time step is chosen to reduce the computational cost of MPC for the chosen prediction horizon. In PI- SINDYc, the parameters for x1 and x3 are estimated well as these only depend on x1, x2 and x3. By contrast, x2 has a larger error in the estimated parameters and consists of erroneous parameters to compensate for the missing information. Different selections of variables have been tested, which generally resulted in poor models, except for the selected combination. The resulting models are generally not sparse, except where a direct relationship exists between variables. This suggests that SINDY indicates direct causal relationships, which can be measured in terms of the sparsity.
Prediction accuracy based on data differing from the training set, but with a similar type of actuation signal, and control results are displayed in figure 12. Both start from an early infection given by x0 = (λ/d, 0.1, 0.1, 0.1, 0.1)T. While a SINDYc model can be identified with near-perfect prediction accuracy, all other models display an error several orders of magnitude larger (see figure 12a). In particular, linear DMDc-based models diverge significantly from the true trajectory for some variables, while capturing the right trend in other variables. The NN and the PI-SINDYc model based on partial state information generally stay closer to, and even temporarily match, the true trajectory. Interestingly, while MPC using PI-SINDYc successfully drives the system to the desired steady-state behaviour, with a slightly larger cost than SINDYc, the NN controller is unable to establish the successful immune response by applying constant treatment (figure 12b). Note that the actuation depends strongly on the prediction and control horizon chosen for the
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

Table 2. Capabilities and challenges of DMDc, SINDYc and NN models. The model with the strongest performance is underlined.
20
property DMDc
training with limited data strong
very few samples are
SINDYc
strong
well suited for low and
NN
weak
requires long time series to
medium amount of data high-dimensionality strong fair strong
can handle high-dim. data limited by the library size in combination with
SVD
………………………………………………………………………………………………………………………………………………………………………………….
sufficient ………………………………………………………………………………………………………………………………………………………………………………….
learn predictive models
nonlinearities weak/fair
linear and weakly
nonlinear, however
strong
suitable for strongly nonlinear systems
strong
suitable for strongly nonlinear systems
with performance loss ………………………………………………………………………………………………………………………………………………………………………………….
prediction performance fair strong strong ………………………………………………………………………………………………………………………………………………………………………………….
Control performance fair strong strong ………………………………………………………………………………………………………………………………………………………………………………….
noise robustness weak strong fair
high sensitivity w.r.t. noise intrinsic robustness due to can handle low noise levels
sparse regression ………………………………………………………………………………………………………………………………………………………………………………….
parameter robustness strong strong weak
high sensitivity w.r.t. initial
weights of the network ………………………………………………………………………………………………………………………………………………………………………………….
training time strong strong weak ………………………………………………………………………………………………………………………………………………………………………………….
execution time strong strong weak fast optimization routines
exist for linear systems ………………………………………………………………………………………………………………………………………………………………………………….
optimization; further analysis has shown that a smaller horizon for the NN controller yields a time-varying, however still unsuccessful, treatment. We varied the number of hidden layers (up to 3), the number of neurons (up to 100), the type of activation function, the number of delays (up to 100) in the state and input variables and the amount of training data (≈ 600 different initial conditions). However, these did not significantly change the performance of the model. The type of data (not just the amount) is particularly critical for training a NN. Designing experiments, i.e. a good forcing signal that explores the system behaviour and yields dynamically rich training data, is a challenge of its own. The linear DMDc and eDMDc models fail too. While the eDMDc model starts with the correct frequency, detrimental treatment is administered thereafter close to the desired state, which gives rise to new growth of infected cells, x2. Interestingly, augmenting the state vector with delay coordinates results in a successful treatment (with performance close to the SINDYc models), in contrast to the strategy to augment the state with nonlinear measurements of the state as in eDMDc.
All models but the NN, which has been trained on a significantly larger amount of data, have been trained on the same amount of data, a single trajectory starting from an initial condition which is relatively far from the desired behaviour. Thus, these models are required to generalize well, i.e. perform well far from the region in which they have been initially trained. Using more data would certainly help to improve the prediction accuracy of some of these models, in particular, if these require a large number of parameters to be estimated. However, this would pose additional challenges in real-time applications with abrupt system changes, as this requires robust model formation and adaptation from few measurements.
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

7. Discussion and conclusion
In conclusion, we have demonstrated the effective integration of data-driven sparse model discovery for MPC in the low-data limit. The sparse identification of nonlinear dynamics (SINDY) algorithm has been extended to discover nonlinear models with actuation and control, resulting in interpretable and parsimonious models. Moreover, because SINDY only identifies the few active terms in the dynamics, it requires less data than many other leading machine learning techniques, such as NNs, and prevents overfitting. When integrated with MPC, SINDY provides computationally tractable and accurate models that can be trained on very little data. The resulting SINDY-MPC framework is capable of controlling strongly nonlinear systems, purely from measurement data, and the model identification is fast enough to discover models in real- time, even in response to abrupt changes to the model. The SINDY-MPC approach is compared with MPC based on data-driven linear models and NN models on four nonlinear dynamical systems of different complexities and challenges: the weakly nonlinear Lotka–Volterra system, the chaotic Lorenz system, the non-affine F8 crusador model, and the HIV/immune response system, which variables are of different order of magnitudes and where only partial state information is available.
The relative strengths and weaknesses of each method are summarized in table 2. By nearly every metric, linear DMDc models and nonlinear SINDYc models outperform NN models. In fact, DMDc may be seen as the limit of SINDYc when the library of candidate terms is restricted to linear terms. SINDY-MPC provides the highest performance control and requires significantly less training data and execution time compared with NN. However, for very low amounts of training data, DMDc provides a useful model until the SINDYc algorithm has enough data to characterize the dynamics. Thus, we advocate the SINDY-MPC framework for effective and efficient nonlinear control, with DMDc as a stopgap after abrupt changes until a new SINDYc model can be identified. Note that a crucial step in SINDY is the choice of library functions, which is often informed by expert knowledge about what category of nonlinearities to include. A poor choice of the library will generally yield a non-sparse model. Without any prior knowledge about the system type, a sweep through different classes of candidate functions is required. However, once a model is learned from a sufficiently rich library, the model is often able to generalize beyond the training data. If the model structure is not fixed, but varies heterogeneously in state space, NNs may provide a more flexible and generalizable architecture to represent the dynamics. A heterogeneous model structure can potentially be incorporated into SINDy by additionally learning a library of models [72,73].
This work motivates a number of future extensions and investigations. Although the preliminary application of SINDYc for MPC is encouraging, this study does not leverage many of the powerful new techniques in sparse model identification. Figure 3 provides a schematic of the modularity and demonstrated extensions that are possible within the SINDy framework. In realistic applications, the system may be extremely high-dimensional, and the SINDy library does not scale well with the size of the data. Fortunately, many high-dimensional systems evolve on a low-dimensional attractor, and it is often possible to identify a model on this attractor, for example by identifying a SINDy model on low-dimensional coordinates obtained through a singular value decomposition [2] or manifold learning [74]. In other applications, full-state measurements are unavailable, and the system must be characterized by limited measurements. It has recently been shown that delay coordinates provide a useful embedding to identify simple models of chaotic systems [53], building on the celebrated Takens embedding theorem [75]. Delay coordinates also define intrinsic coordinates for the Koopman operator [53], which provides a simple linear embedding of nonlinear systems [76,77]. Koopman models have recently been used for MPC [24,25] and have been identified using SINDy regression [78] and subsequently used for optimal control [78]. Recently, SINDY has been extended to modify an existing model based on new incoming measurements to enable rapid model recovery from abrupt changes to the system [30]. Learning quickly from limited measurements is an important task, which may be viewed in terms of design of experiments; specifically, optimizing the actuation input to collect
21
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

the most informative measurements to learn a more predictive model faster. This would require the formulation of a different cost function, which measures the predictive power of the model, to determine future actuation inputs. Rapid learning is also related to the question of quantity versus quality of data and identifiability [48,49]; more data is usually better, although it is possible to work with less data if it is representative of the system. Further, similar methods could be used to optimize sensors and exploit partial measurements within the SINDY-MPC framework. All of these innovations suggest a shift from the perspective of big data to the control-oriented perspective of smart data.
Figure 3 also demonstrates innovations to the SINDy regression to include physical constraints, known model structure, and model selection, which may all benefit the goal of real- time identification and control. Known symmetries, conservation laws, and constraints may be readily included in both the SINDYc and DMDc modelling frameworks [31], as they are both based on least-squares regression, possibly with sequential thresholding. It is thus possible to use a constrained least-squares algorithm, for example, to enforce energy conserving constraints in a fluid system, which manifest as anti-symmetric quadratic terms [31]. Enforcing constraints has the potential to further reduce the amount of data required to identify models, as there are less free parameters to estimate, and the resulting systems have been shown to have improved stability in some cases. It is also possible to extend the SINDy algorithm to identify models in libraries that encode richer dynamics, such as rational function nonlinearities [79]. Finally, incorporating information criteria provides an objective metric for model selection among various candidate SINDy models with a range of complexity.
The SINDY-MPC framework has significant potential for the real-time control of strongly nonlinear systems. Moreover, the rapid training and execution times indicate that SINDy models may be useful for rapid model identification in response to abrupt model changes, and this warrants further investigation. The ability to identify accurate and efficient models with small amounts of training data may be a key enabler of recovery in time-critical scenarios, such as model changes that lead to instability. In addition, for broad applicability and adoption, the SINDy modelling framework must be further investigated to characterize the effect of noise, derive error estimates, and provide conditions and guarantees of convergence. These future theoretical and analytical extensions are necessary to certify the model-based control performance.
Data accessibility. The code used in this work is made available at: https://github.com/eurika-kaiser/ SINDY-MPC. The data can be generated using the code.
Authors’ contributions. All authors conceived of the work, designed the study and drafted the manuscript. E.K. carried out the computations.
Competing interests. We declare we have no competing interests.
Funding. E.K. gratefully acknowledges support by the Washington Research Foundation, the Gordon and Betty Moore Foundation (Award no. 2013-10-29), the Alfred P. Sloan Foundation (Award no. 3835), and the University of Washington eScience Institute. S.L.B. and J.N.K. acknowledge support from the Defense Advanced Research Projects Agency (DARPA contract HR011-16-C-0016 and PA-18-01-FP-125). S.L.B. acknowledges support from the Army Research Office (W911NF-17-1-0306 and W911NF-17-1-0422). J.N.K. acknowledges support from the Air Force Office of Scientific Research (FA9550-17-1-0329).
Acknowledgements. The authors gratefully acknowledge many valuable discussions with Josh Proctor.
References
1. Kutz JN, Brunton SL, Brunton BW, Proctor JL. 2016 Dynamic mode decomposition: data-driven modeling of complex systems. Philadelphia, PA: SIAM.
2. Brunton SL, Proctor JL, Kutz JN. 2016 Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc. Natl Acad. Sci. 113, 3932–3937. (doi:10.1073/pnas.1517384113)
3. Allgöwer F, Badgwell TA, Qin JS, Rawlings JB, Wright SJ. 1999 Nonlinear predictive control and moving horizon estimation: an introductory overview. In Advances in control: Highlights of ECC’99 (ed. PM Frank), pp. 391–449. London, UK: Springer.
22
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

4. Camacho EF, Alba CB. 2013 Model predictive control. Berlin, Germany: Springer Science & Business Media.
5. Skogestad S, Postlethwaite I. 2005 Multivariable feedback control: analysis and design, 2nd edn. Hoboken, NJ: John Wiley & Sons, Inc.
6. Dullerud GE, Paganini F. 2000 A course in robust control theory: a convex approach. Texts in Applied Mathematics. Berlin, Germany: Springer.
7. Garcia CE, Prett DM, Morari M. 1989 Model predictive control: theory and practice — a survey. Automatica 25, 335–348. (doi:10.1016/0005-1098(89)90002-2)
8. Morari M, Lee JH. 1999 Model predictive control: past, present and future. Comput. Chem. Eng. 23, 667–682. (doi:10.1016/S0098-1354(98)00301-9)
9. Lee JH. 2011 Model predictive control: review of the three decades of development. Int. J. Control Autom. Syst. 9, 415–424. (doi:10.1007/s12555-011-0300-6)
10. Mayne DQ. 2014 Model predictive control: recent developments and future promise. Automatica 50, 2967–2986. (doi:10.1016/j.automatica.2014.10.128)
11. Eren U, Prach A, Koçer BB Rakovic ́ SV, Kayacan B, Açıkmes ̧e E. 2017 Model predictive control in aerospace systems: current state and opportunities. J. Guid. Control Dyn. 40, 1541–1566. (doi:10.2514/1.G002507)
12. Brunton SL, Noack BR. 2015 Closed-loop turbulence control: progress and challenges. Appl. Mech. Rev. 67, 050801–1–050801–48. (doi:10.1115/1.4031175)
13. Juang JN, Pappa RS. 1985 An eigensystem realization algorithm for modal parameter identification and model reduction. J. Guid. Control Dyn. 8, 620–627. (doi:10.2514/3.20031)
14. Brockett RW. 1976 Volterra series and geometric control theory. Automatica 12, 167–176.
(doi:10.1016/0005-1098(76)90080-7)
15. Boyd S, Chua LO, Desoer CA. 1984 Analytical foundations of Volterra series. IMA J. Math.
Control Inf. 1, 243–282. (doi:10.1093/imamci/1.3.243)
16. Maner BR, Doyle FJ, Ogunnaike BA, Pearson RK. 1994 A nonlinear model predictive control
scheme using second order Volterra models. In American Control Conference, Baltimore, MD, 29
June–1 July, vol. 3, pp. 3253–3257. Piscataway, NJ: IEEE.
17. Akaike H. 1969 Fitting autoregressive models for prediction. Ann. Inst. Stat. Math. 21, 243–247.
(doi:10.1007/BF02532251)
18. Billings SA. 2013 Nonlinear system identification: NARMAX methods in the time, frequency, and
spatio-temporal domains. Chichester, UK: John Wiley & Sons.
19. Lippmann R. 1987 An introduction to computing with neural nets. IEEE Assp Mag. 4, 4–22.
(doi:10.1109/MASSP.1987.1165593)
20. Draeger A, Engell S, Ranke H. 1995 Model predictive control using neural networks. IEEE
Control Syst. Mag. 15, 61–66. (doi:10.1109/37.466261)
21. Wang T, Gao H, Qiu J. 2016 A combined adaptive neural network and nonlinear model
predictive control for multirate networked industrial process control. IEEE Trans. Neural Netw.
Learn. Syst. 27, 416–425. (doi:10.1109/TNNLS.2015.2411671)
22. Aggelogiannaki E, Sarimveis H. 2008 Nonlinear model predictive control for distributed
parameter systems using data driven artificial neural network models. Comput. Chem. Eng.
32, 1225–1237. (doi:10.1016/j.compchemeng.2007.05.002)
23. Williams MO, Kevrekidis IG, Rowley CW. 2015 A data-driven approximation of the
Koopman operator: extending dynamic mode decomposition. J. Nonlin. Sci. 25, 1307–1346.
(doi:10.1007/s00332-015-9258-5)
24. Korda M, Mezic ́ I. 2016 Linear predictors for nonlinear dynamical systems: Koopman operator
meets model predictive control. (http://arxiv.org/abs/1611.03537)
25. Peitz S, Schäfer K, Ober-Blöbaum S, Eckstein J, Köhler U, Dellnitz M. 2016 A multiobjective
MPC approach for autonomously driven electric vehicles. (http://arxiv.org/abs/1610.08777)
26. Peng H, Wu J, Inoussa G, Deng Q, Nakano K. 2009 Nonlinear system modeling and predictive control using the RBF nets-based quasi-linear ARX model. Control Eng. Pract. 17, 59–66.
(doi:10.1016/j.conengprac.2008.05.005)
27. Zhang T, Kahn G, Levine S, Abbeel P. 2016 Learning deep control policies for autonomous
aerial vehicles with MPC-guided policy search. In IEEE International Conference on Robotics and
Automation, Stockholm, Sweden, 16–21 May, pp. 528–535. Piscataway, NJ: IEEE.
28. Bongard J, Lipson H. 2007 Automated reverse engineering of nonlinear dynamical systems.
Proc. Natl Acad. Sci. USA 104, 9943–9948. (doi:10.1073/pnas.0609476104)
29. Schmidt M, Lipson H. 2009 Distilling free-form natural laws from experimental data. Science
23
324, 81–85. (doi:10.1126/science.1165893)
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

30. Quade M, Abel M, Kutz JN, Brunton SL. 2018 Sparse identification of nonlinear dynamics for rapid model recovery. Chaos 28, 063116. (doi:10.1063/1.5027470)
31. Loiseau JC, Brunton SL. 2018 Constrained sparse Galerkin regression. J. Fluid Mech. 838, 42–67. (doi:10.1017/jfm.2017.823)
32. Proctor JL, Brunton SL, Kutz JN. 2016 Dynamic mode decomposition with control. SIAM J. Appl. Dyn. Syst. 15, 142–161. (doi:10.1137/15M1013857)
33. Li K, Peng JX, Irwin GW. 2005 A fast nonlinear model identification method. IEEE Trans. Autom. Control 50, 1211–1216. (doi:10.1109/TAC.2005.852557)
34. Chen T, Andersen MS, Ljung L, Chiuso A, Pillonetto G. 2014 System identification via sparse multiple kernel-based regularization using sequential convex optimization techniques. IEEE Trans. Autom. Control 59, 2933–2945. (doi:10.1109/TAC.2014.2351851)
35. Xu W, Bai EW, Cho M. 2014 System identification in the presence of outliers and random noises: a compressed sensing approach. Automatica 50, 2905–2911. (doi:10.1016/j.automatica. 2014.10.017)
36. Pan W, Yuan Y, Gonçalves J, Stan GB. 2012 Reconstruction of arbitrary biochemical reaction networks: a compressive sensing approach. In 2012 IEEE 51st Annual Conference on Decision and Control (CDC), Maui, HI, 10–13 December, pp. 2334–2339. Piscataway, NJ: IEEE.
37. Calafiore GC, El Ghaoui LM, Novara C. 2015 Sparse identification of posynomial models. Automatica 59, 27–34. (doi:10.1016/j.automatica.2015.06.003)
38. Pan W, Yuan Y, Gonçalves J, Stan GB. 2016 A sparse Bayesian approach to the identification of nonlinear state-space systems. IEEE Trans. Autom. Control 61, 182–187. (doi:10.1109/TAC.2015.2426291)
39. Nelles O. 2013 Nonlinear system identification: from classical approaches to neural networks and fuzzy models. Berlin, Germany: Springer Science & Business Media.
40. Pillonetto G, Dinuzzo F, Chen T, De Nicolao G, Ljung L. 2014 Kernel methods in system identification, machine learning and function estimation: a survey. Automatica 50, 657–682. (doi:10.1016/j.automatica.2014.01.001)
41. Rudin LI, Osher S, Fatemi E. 1992 Nonlinear total variation based noise removal algorithms. Physica D 60, 259–268. (doi:10.1016/0167-2789(92)90242-F)
42. Chartrand R. 2011 Numerical differentiation of noisy, nonsmooth data. ISRN Appl. Math. 2011, 1–11. (doi:10.5402/2011/164564)
43. Tibshirani R. 1996 Regression shrinkage and selection via the lasso. J. R. Stat. Soc. B 58, 267–288.
44. Tropp JA. 2006 Just relax: Convex programming methods for identifying sparse signals in noise. IEEE Trans. Inf. Theory 52, 1030–1051. (doi:10.1109/TIT.2005.864420)
45. Su W, Bogdan M, Candès EJ. 2016 False discoveries occur early on the Lasso path. (http://arxiv.org/abs/1511.01957)
46. Zhang L, Schaeffer H. 2018 On the convergence of the SINDy algorithm. (http://arxiv.org/ abs/1805.06445)
47. Zheng P, Askham T, Brunton SL, Kutz JN, Aravkin AY. 2018 A unified framework for sparse relaxed regularized regression: SR3. (http://arxiv.org/abs/1807.05411)
48. Gevers M, Bazanella AS, Coutinho DF, Dasgupta S. 2013 Identifiability and excitation of polynomial systems. In 52nd IEEE Conference Decision and Control (CDC), Florence, Italy, 10–13 December, pp. 4278–4283. IEEE.
49. Alkhoury Z, Petreczky M, Mercère G. 2017 Identifiability of affine linear parameter-varying models. Automatica 80, 62–74. (doi:10.1016/j.automatica.2017.01.029)
50. Mangan NM, Kutz JN, Brunton SL, Proctor JL. 2017 Model selection for dynamical systems via sparse regression and information criteria. Proc. R. Soc. A 473, 1–16. (doi:10.1098/ rspa.2017.0009)
51. Rudy SH, Brunton SL, Proctor JL, Kutz JN. 2017 Data-driven discovery of partial differential equations. Sci. Adv. 3, e1602614. (doi:10.1126/sciadv.1602614)
52. Schaeffer H. 2017 Learning partial differential equations via data discovery and sparse optimization. Proc. R. Soc. A 473, 20160446. (doi:10.1098/rspa.2016.0446)
53. Brunton SL, Brunton BW, Proctor JL, Kaiser E, Kutz JN. 2017 Chaos as an intermittently forced linear system. Nat. Commun. 8, 1–9. (doi:10.1038/s41467-017-00030-8)
54. Tran G, Ward R. 2016 Exact recovery of chaotic systems from highly corrupted data. (http://arxiv.org/abs/1607.01067)
55. Schaeffer H, McCalla SG. 2017 Sparse model selection via integral terms. Phys. Rev. E 96, 023302. (doi:10.1103/PhysRevE.96.023302)
24
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

56. Rowley CW, Mezic ́ I, Bagheri S, Schlatter P, Henningson D. 2009 Spectral analysis of nonlinear flows. J. Fluid Mech. 645, 115–127. (doi:10.1017/S0022112009992059)
57. Schmid PJ. 2010 Dynamic mode decomposition of numerical and experimental data. J. Fluid Mech. 656, 5–28. (doi:10.1017/s0022112010001217)
58. Tu JH, Rowley CW, Luchtenburg DM, Brunton SL, Kutz JN. 2014 On dynamic mode decomposition: theory and applications. J. Comput. Dyn. 1, 391–421. (doi:10.3934/ jcd.2014.1.391)
59. Holmes PJ, Lumley JL, Berkooz G, Rowley CW. 2012 Turbulence, coherent structures, dynamical systems and symmetry, 2nd edn. Cambridge, England: Cambridge University Press.
60. Johansen TA. 2011 Introduction to nonlinear model predictive control and moving horizon estimation. In Selected topics on constrained and nonlinear control, pp. 187–240. Bratislava, Slovakia: STU and Trondheim, Norway: NTNU.
61. Lee S, Lee D, Oh H. 2005 Technological forecasting at the Korean stock market: a dynamic competition analysis using Lotka-Volterra model. Technol. Forecast. Soc. Change 72, 1044–1057. (doi:10.1016/j.techfore.2002.11.001)
62. Venturino E. 1994 The influence of diseases on Lotka-Volterra systems. Rocky Mt. J. Math. 24, 381–402. (doi:10.1216/rmjm/1181072471)
63. Hornik K, Stinchcombe M, White H. 1989 Multilayer feedforward networks are universal approximators. Neural Netw. 2, 359–366. (doi:10.1016/0893-6080(89)90020-8)
64. Schroeder M. 1970 Synthesis of low-peak-factor signals and binary sequences with low autocorrelation (Corresp.). IEEE Trans. Inf. Theory 16, 85–89. (doi:10.1109/TIT.1970.1054411)
65.Lorenz EN. 1963 Deterministic nonperiodic flow. J. Atmos. Sci. 20, 130–141. (doi:10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2)
66. Garrard WL, Jordan JM. 1977 Design of nonlinear automatic flight control systems. Automatica 13, 497–505. (doi:10.1016/0005-1098(77)90070-X)
67. Çimen T, Banks SP. 2004 Global optimal feedback control for general nonlinear systems with nonquadratic performance criteria. Syst. Control Lett. 53, 327–346. (doi:10.1016/j.sysconle. 2004.05.008)
68. Yan Z, Wang J. 2012 Model predictive control of nonlinear systems with unmodeled dynamics based on feedforward and recurrent neural networks. IEEE Trans. Ind. Inform. 8, 746–756. (doi:10.1109/TII.2012.2205582)
69. Perelson AS, Nelson PW. 1999 Mathematical analysis of HIV-1 dynamics in vivo. SIAM Rev. 41, 3–44. (doi:10.1137/S0036144598335107)
70. Zurakowski R, Teel AR. 2006 A model predictive control based scheduling method for HIV therapy. J. Theor. Biol. 238, 368–382. (doi:10.1016/j.jtbi.2005.05.004)
71. Wodarz D. 2001 Helper-dependent vs. helper-independent CTL responses in HIV infection: implications for drug therapy and resistance. J. Theor. Biol. 213, 447–459. (doi:10.1006/jtbi.2001.2426)
72. Brunton SL, Tu JH, Bright I, Kutz JN. 2014 Compressive sensing and low-rank libraries for classification of bifurcation regimes in nonlinear dynamical systems. SIAM J. Appl. Dyn. Syst. 13, 1716–1732. (doi:10.1137/130949282)
73. Sargsyan S, Brunton SL, Kutz JN. 2015 Nonlinear model reduction for dynamical systems using sparse sensor locations from learned libraries. Phys. Rev. E 92, 033304. (doi:10.1103/PhysRevE.92.033304)
74. Loiseau JC, Noack BR, Brunton SL. 2018 Sparse reduced-order modelling: sensor-based dynamics to full-state estimation. J. Fluid Mech. 844, 459–490. (doi:10.1017/jfm.2018.147)
75. Takens F. 1981 Detecting strange attractors in turbulence. Lect. Notes Math. 898, 366–381. (doi:10.1007/BFb0091924)
76. Mezic ́ I. 2005 Spectral properties of dynamical systems, model reduction and decompositions. Nonlinear Dyn. 41, 309–325. (doi:10.1007/s11071-005-2824-x)
77. Budišic ́ M, Mohr R, Mezic ́ I. 2012 Applied Koopmanism. Chaos 22, 047510. (doi:10.1063/ 1.4772195)
78. Kaiser E, Kutz JN, Brunton SL. 2017 Data-driven discovery of Koopman eigenfunctions for control. (http://arxiv.org/abs/1707.01146).
79. Mangan NM, Brunton SL, Proctor JL, Kutz JN. 2016 Inferring biological networks by sparse identification of nonlinear dynamics. IEEE Trans. Mol. Bio. Multi-Scale Commun. 2, 52–63. (doi:10.1109/TMBMC.2016.2633265)
25
rspa.royalsocietypublishing.org Proc. R. Soc. A 474: 20180335 ……………………………………………

Related Posts