NAMICAL SYSTEMS USING DEEP NEURAL NETWORK
Year: 2021

Abstract:
A novel nonlinear system identification method that produces state estimator and predictor directly usable for model predictive control (MPC) is proposed in this paper. The main feature of the proposed method is that it uses a neural network with a bottleneck layer between the state estimator and predictor to represent the input-output dynamics, and it is proven that the state of the dynamical system can be extracted from the bottleneck layer based on the observability of the target system. The training of the network is shown to be a natural nonlinear extension of the subspace state-space system identification method established for linear dynamical systems. This correspondence gives interpretability to the resulting model based on linear control theory. The usefulness of the proposed method and the interpretability of the model are demonstrated through an illustrative example of MPC.

heading: 1 INTRODUCTION
text: Model predictive control (MPC) has been widely used as a practical method to control nonlinear dynamical systems (Qin & Badgwell, 2000; Vazquez et al., 2014; Karamanakos et al., 2020). As the name implies, MPC requires a prediction model, and in most cases, state-space models are used. When we want to apply MPC to a complex system whose structure and state are unknown, model construction from input-output data, i.e., system identification, is required. For linear state-space system identification, the subspace methods (Favoreel et al., 2000; Katayama, 2005; van der Veen et al., 2013) are widely used. Its fundamental idea is to reduce the dimensionality of the measured past and future input and output data, and it strongly depends on the linearity of the system. Although there have been continuous attempts to extend the subspace methods to nonlinear systems with specific structures such as Wiener and Hammerstein structures (e.g., Verhaegen & Westwick, 1996; Goethals et al., 2005; Katayama & Ase, 2016) and related research has also been conducted in the field of machine learning as spectral learning of Hidden Markov Models (Kawahara et al., 2007; Hsu et al., 2012; Kandasamy et al., 2016), there is still a lot of room for research in extending the subspace method to general nonlinear dynamical systems.
On the other hand, recent advances in machine learning techniques have had an impact on the research field of system identification (Ljung et al., 2020a;b). For example, deep state space models (Rangapuram et al., 2018), and Kalman variational auto-encoder (Fraccaro et al., 2017) apply deep learning to state space models. Among those technologies, autoencoder is a notable dimensional reduction or feature extraction technique with deep neural networks, whose traits derives from topological constraints of own bottleneck structure. Although several pieces of research apply it to state-space system identification (Masti & Bemporad, 2018; Beintema et al., 2021), these works do not provide a theoretical basis for obtaining the state sequence.
This paper aims to propose a theoretically supported and interpretable nonlinear dynamical system identification method by simply extending the subspace method with neural networks. For the proposed method, it can be shown based on the observability of the target system that when a neural network with an appropriate bottleneck is trained, the state estimator and predictor are formed
across the bottleneck layer. The obtained model can be interpreted based on the knowledge from linear system identification methods, and this aspect is explained through an illustrative example.
In this paper, the set of integers and the natural numbers are denoted by Z and N, respectively. The set {sa, sa+1, . . . , sb−1, sb} is denoted as {sk}bk=a for conciseness. The n ×m zero matrix is denoted as 0n×m, and 0n := 0n×1. For a matrix M , M † is the Moore–Penrose inverse of M .

heading: 2 PROBLEM FORMULATION & NOTATION
text: We consider the following discrete-time nonlinear dynamical system:
xt+1 = f(xt,ut), yt = h(xt) (1)
with t ∈ Z the time index, xt ∈ Rnx the state, yt ∈ Rny the output, ut ∈ Rnu the input, f : Rnx × Rnu → Rnx the state-transition map and h : Rnx → Rny the output map. As in the usual setting for system identification, xt is a latent variable and ut and yt are measurable with noise. In the following discussion, the fundamental properties of the proposed method for nonlinear target systems are presented in the absence of measurement noise and disturbances. The effects of noise and disturbance are indirectly provided for linear target systems in section 4 by showing that the proposed method is equivalent to the subspace identification method.
For convenience, We define the sequence of input and output signals in the past and future from a certain time t as follows:
upt :=  ut−hp
ut−hp+1 ...
ut−1
 , ypt :=  yt−hp yt−hp+1
... yt−1
 , uft :=  ut ut+1
... ut+hf−1
 , yft :=  yt yt+1
... yt+hf−1  (2) where the superscripts p and f indicate past and future, and hp ∈ N and hf ∈ N are horizons for corresponding directions, respectively. Using these notations, we formulate the problem as follows.
Problem 1. Given the measured input-output data {(ut,yt)}T+hf−1t=−hp+1, and the design parameters hp,hf , nx̂ ∈ N, construct a model that consists of a state estimator E, which maps (upt ,y p t ) 7→ x̂t, and a predictor P , which maps (x̂t,uft) 7→ yft . Here, x̂t ∈ Rnx̂ is a state equivalent to xt but in an arbitrary coordinate system and T ∈ N indicates the size of the dataset. Remark 1. In the usual state-space system identification, the goal is to construct a model of f and h, i.e., a state-space model. However, in many applications, including model predictive control, the state estimator and multi-step ahead predictor are often reconstructed from the resulting state-space model. In particular, for nonlinear systems that are difficult to analyze, the merits of going through state-space models are questionable. Therefore, we focus here on the construction of the state estimator E and predictor P and only briefly discuss the identification of f and h in Appendix A.
To guarantee Problem 1 to be solvable, we make some assumptions. With the vectors defined above (2), we can write the system (1) as
xt = f hp ( xt−hp ,u p t ) , ypt = h hp ( xt−hp ,u p t ) (3)
where
fk ( xt, [ u>t , . . . ,u > t+k ]>) := f (· · · f (f (xt,ut) ,ut+1) · · · ,ut+k) (4a)
hk ( xt, [ u>t , . . . ,u > t+k ]>) :=  h(xt) h ◦ f(xt,ut) ...
h ◦ fk−1 ( xt, [ u>t , . . . ,u > t+k−1
]>)  . (4b)
And define the uniform k-observability as follows. Definition 1 (uniform k-observability (Moraal & Grizzle, 1995)). If the mapping
Rnx × (Rnu)k → (Rny )k × (Rnu)k by (x,u) 7→ ( hk(x,u),u ) (5)
is injective, the system (1) is said to be uniformly k-observable (Moraal & Grizzle, 1995).
Then, we make the following assumptions and theorem.
Assumption 1. The system (1) is uniformly k-observable, where k = min (hp,hf). Remark 2. According to Stark et al. (1997), typical dynamical systems are k-observable for k ≥ 2nx + 1.
Assumption 2. The dimension of the state estimate is large enough, that is, nx̂ ≥ nx. Theorem 1. Under Assumptions 1 and 2, Problem 1 has a solution.
Proof. From Assumption 1, there exists the map Φhp : {( up, hhp(xp,up) ) ∣∣ xp ∈ Rnx ,up ∈ (Rnu)hp}→ Rnx (6)
that maps (upt ,y p t ) 7→ xt−hp , and the map
Ψ : {( up, hhp(xp,up) ) ∣∣ xp ∈ Rnx ,up ∈ (Rnu)hp}→ Rnx (7)
that maps (upt ,y p t ) 7→ fhp ( Φhp (upt ,y p t ) ,u p t ) also exists. Since the pair of maps Ψ and hhf is a pair of state estimator and predictor that is a solution to Problem 1, it is guaranteed that there is a solution when nx̂ ≥ nx.

heading: 3 PROPOSED METHOD
text: This section presents the proposed nonlinear system identification approach. Also, nonlinear MPC based on the obtained model is described.

heading: 3.1 TRAINING STATE ESTIMATOR & PREDICTOR
text: In the proposed method, we train the neural network which has the structure of figure 1, where the state estimator E : (Rnu)hp × (Rny )hp → Rnx̂ maps the past input and output to the current state and the predictor P : Rnx̂ × Rhfnu → Rhf×ny maps the current state and the future input to the future output.
We set nx̂ < hpny for the bottleneck structure of the network to take a lower-dimensional state, and the training problem is formulated as the following minimization problem
min E,P
1
T T∑ t=1 ∥∥∥yft − P (E (upt ,ypt )︸ ︷︷ ︸ estimate of xt ,uft )
︸ ︷︷ ︸ prediction of yft
∥∥∥2. (8)
Since the state xt cannot be measured, the state estimator E and predictor P are coupled in series and trained in the same way as the encoder and decoder in autoencoder. To solve (8), we can utilize modern batch optimization algorithms for deep learning (e.g. Adam by Kingma & Ba, 2014).
Here, we show that x̂t := E (u p t ,y p t ) is also a state of the system in the ideal situation.
Theorem 2. Assume the training (8) is perfectly done, that is, for all xp ∈ Rnx , up ∈ (Rnu)hp , uf ∈ (Rnu)hf , the equation
P ( E (up,yp) ,uf ) = yf , (9)
where yf ( xp,up,uf ) := hhf ( x,uf ) , x (xp,up) := fhp (xp,up), yp (xp,up) := hhp (xp,up) and arguments are omitted for conciseness, holds. Then, for all xp ∈ Rnx and up ∈ (Rnu)hp , E ( up, hhp (xp,up) ) is a state equivalent to fhp (xp,up).
Proof. From hf -observability of the system, the map Φhf : {( uf , hhf ( x,uf )) ∣∣ x ∈ Rnx ,uf ∈ (Rnu)hf}→ Rnx that maps (uf ,yf) 7→ x (10) exists, and the map
Γ : Rnx̂ × (Rnu)hf → Rnx that maps ( x̂,uf ) 7→ Φhf ( uf , P ( x̂,uf )) (11)
also exists. Then, a map obtained by fixing the second argument of Γ, e.g., Γ0(x̂) := Γ (x̂,0hfnu) (12) maps x̂ to x, and x̂ is proved to be a state.
Remark 3. The existence of the map from x̂ to x is ensured by the predictor P trained with a common input 0hfnu for all states as in (12). Although this is just one way of constructing a map that is possible under the strong assumption of Theorem 2, this implies that to obtain an appropriate state estimator, the training data must be rich enough to include situations where common inputs are applied in diverse states. Corollary 1. When the assumption of Theorem 2 holds, P is the predictor that predicts the sequence of future output from a state and a sequence of future input.

heading: 3.2 MODEL PREDICTIVE CONTROL WITH ESTIMATOR & PREDICTOR
text: Next, the use of the obtained model in MPC is briefly described.
Here we reuse the symbols for the past and future input/output sequences upt ,y p t ,u f t,y f t defined in (2), but note that the MPC is done after the training and these input-output signals are not the ones at the time of training.
In MPC, a loss function lt : Rnu × Rny → R that evaluates the input and output at each time t is given. And, at each time step t, the optimization problem
minimize uf∈(Rnu )hf hf∑ k=1 lt+k−1 (( uf ) (k−1)nu+1:knu , ( yf ) (k−1)ny+1:kny ) (13)
s.t. yf = P ( E (upt ,y p t ) ,u f )(
uf ,yf ) satisfies constraints on the control system
is solved to obtain uf , which is a candidate for the future input sequence uft. Then, the candidate for most immediate future input [ Inu ,0nu×nu(hf−1) ] uf is actually applied to the system and moves to the next step t+ 1.
Usually, state-space models are used in MPC, and a state estimator like extended Kalman filter and a predictor are constructed from the models. On the other hand, the state estimator and predictor obtained by the proposed method can be directly used for MPC. The state estimator obtained by the proposed method is a static mapping from the previous inputs and outputs on the finite horizon hp, and there is no risk of divergence, unlike commonly used filters with infinite impulse response (IIR) structures. This is an important advantage when using the state estimator for MPC.

heading: 4 INTERPRETATION OF PROPOSED METHOD AND MODEL
text: The proposed method can be viewed as a natural extension of the subspace methods for linear systems, and the similarity with the subspace methods gives interpretability to the model obtained by the proposed method. In this section, we briefly explain the framework of the subspace method and describe how the model obtained by the proposed method can be interpreted based on it.

heading: 4.1 RELATION TO SUBSPACE METHODS FOR LINEAR DYNAMICAL SYSTEMS
text: Here, we consider restricting the state estimator E and predictor P of the proposed method to linear maps and representing them by matrices E ∈ Rnx̂×hp(nu+ny) and P ∈ Rhfny×(nx+hfnu). When the target system is linear, the optimal state estimator and predictor are known to be linear maps, and can be represented by these matrices if hp and hf are sufficiently large. In this case, the training of the state estimator and predictor (8) is reduced to the following nonlinear least-squares problem:
min E∈Rnx̂×hp(nu+ny)
Px∈Rhfny×nx̂ Pu∈Rhfny×hfnu
1
T T∑ t=1 ∥∥∥∥yft − PxE [uptypt ] − Puuft ∥∥∥∥2 , (14) where we split P = [Px,Pu] for convenience. Also, by defining the matrix Up := [u p 1 ,u p 2 , . . . ,u p T ] and the matrices Y p, U f , Y f in the same manner, we can rewrite the objective function as follows
min E∈Rnx̂×hp(nu+ny)
Px∈Rhfny×nx̂ Pu∈Rhfny×hfnu
1
T
∥∥∥∥Y f − PxE [UpY p ] − PuU f ∥∥∥∥2 F . (15)
On the other hand, while the actual computational procedure is more cumbersome due to the use of RQ decomposition for numerical benefits, the outline of the subspace methods can be summarized as follows (Peternell et al., 1996):
(S1) Solve linear least squares problem
min P̂xE∈Rhfny×hp(nu+ny)
Pu∈Rhfny×hfnu
1
T
∥∥∥∥Y f − P̂xE [UpY p ] − PuU f ∥∥∥∥2 F
(16)
to obtain P̂xE ∈ Rhfny×hp(nu+ny) and Pu ∈ Rhfny×hfnu . (S2) Obtain Px and E rank-nx̂ approximation of P̂xE via a ‘weighted’ SVD. (S3) Calculate state-space representation based on E.
As can be seen from the comparison of the two optimization problems (15) and (16), the subspace method is characterized by the fact that instead of solving the nonlinear least-squares problem involving the product of two matrix parameters Px and E, it solves the linear least-squares problem with the product P̂xE as the parameter, and then obtain Px and E via the low-rank approximation in (S2).
The distinct advantage of the procedure (S1)-(S2) is that both linear least-squares and low-rank approximations by SVD yield globally optimal solutions. On the other hand, the solution obtained by (S1)-(S2) is different from the solution of the original optimization problem (15), and there are many variants that have different statistical properties for this misalignment (see e.g., Van Overschee & De Moor, 1994; Peternell et al., 1996; Favoreel et al., 2000).
Therefore, the proposed method can be regarded as a variant of the subspace methods that revert to the underlying nonlinear least-squares problem (15) and extend it to nonlinear systems by replacing linear mappings with neural networks.
Appendix B confirms that the proposed method indeed yields the same results as the ideal subspace identification method for a linear system example. Next, we discuss the interpretation of the obtained model based on this agreement.

heading: 4.2 INTERPRETATION OF MODEL
text: It is also possible to interpret the obtained model based on the knowledge in linear system identification. Only in this section, we assume that the target system is the following linear dynamical system with noise
xt+1 = Axt + But + wt, yt = Cxt + vt, (17)
where A ∈ Rnx×nx , B ∈ Rnx×nu , C ∈ Rny×nx are the system matrices; wt ∼ N (0nu×1,Q) and vt ∼ N ( 0ny×1,R ) are process noise and measurement noise, respectively. When the number of states in the model is accurate (nx̂ = nx) and the number of data T and horizons hp,hf are sufficiently large, the result of steps (S1) and (S2) is consistent with the solution of (15), and what we get as a state estimator E is known to be the Kalman filter for the target system. And, the matrix E becomes the Markov parameters of the Kalman filter, that is,
E = [Eu Ey] , (18) Eu = [ Āhp−1B, . . . , Ā2B, ĀB,B ] , Ey = [ Āhp−1K, . . . , Ā2K, ĀK,K ] , (19)
Ā := A−KC, K := AΣ̂C> ( CΣ̂C> + R )−1 (20)
where K ∈ Rnx×ny is the Kalman gain and Σ̂ is the solution to the discrete-time algebraic Riccati equation
Σ̂ = AΣ̂A> −AΣ̂C> ( CΣ̂C> + R )−1 CΣ̂A> + Q. (21)
From the Kalman filter designed for the target system, we can retrieve the frequency bands that are important for state estimation or the noise characteristics in the target system. Conversely, we can also evaluate whether the obtained state estimator is reasonable or not, based on the knowledge of the target system. For the state estimator E implemented by a neural network, Eu and Ey correspond to ∂E∂up and ∂E ∂yp at the operating point, respectively. And by calculating these, the characteristics of the obtained state estimator can be evaluated.
On the other hand, the predictor Px and Pu are an extended observability matrix and the Markov parameters of the target system, that is,
Px =  C CA
... CAhf−1
 , Pu =  O O · · · · · · O CB O · · · CAB CB . . . ... ... ... . . .
CAhf−2B CAhf−3B · · · CB O
 , O := 0ny×nu . (22)
and characterize the unforced response and input-output characteristics of the target system. For nonlinear predictors, Px and Pu correspond to ∂P∂x̂ and ∂P ∂yp , respectively, from which the characteristics of the target system at the operating point can be illustrated.

heading: 5 ILLUSTRATIVE EXAMPLE
text: Next, we demonstrate the proposed method through a toy example, where the identification and MPC of a cascaded tanks system (Wigren & Schoukens, 2013) is performed.
The state equation of the system is as follows:
xt+1 = [ xt+1,1 xt+1,2 ] = [ max ( xt,1 − k1 √ xt,1 + k2 (ut + wt) , 0 ) max ( xt,2 + k3 √ xt,1 − k4 √ xt,2, 0 ) ] , yt = xt,2 + vt, (23) where nx = 2, nu = ny = 1, k1 = 0.5, k2 = 0.4, k3 = 0.2, k4 = 0.3, wt ∼ N ( 0, 0.352 ) , and
vt ∼ N ( 0, 0.152 ) . The state variables xt,1 and xt,2 are the levels of the upper and the lower tank, and ut denotes the applied voltage to the pump at t. To obtain the training dataset, the system is excited by a sequence of random input
ut = min ( max ( u0t , 0 ) , 5 ) , u0t ∼ N (2, 52), (24)
and a sequence of data {(ut, yt)}1000019t=−19 is obtained. On training the state estimator E and predictor P , We choose the hyperparameters of nx̂ = 2, hp = hf = 20. E and P are implemented by neural networks with two and three hidden layers, respectively, where each layer contains 64 ReLU neurons and every layers are fully connected. For
training, we used the Adam optimizer in PyTorch (Paszke et al., 2019) with the default settings. To prevent overfitting, the 90 % of the samples were used for training, and the training was continued until the fit to the remaining 10 % of the validation data did not improve for 1000 consecutive iterations.
Then, the MPC controller which solves the optimization problem (13) with the cost function
lt(ut, yt) := (yt − rt)2 + 0.05 (ut − ut−1)2 (25) and the constraint 0 ≤ ut ≤ 5 at each time t is implemented with the obtained state estimator and predictor. Here, rt ∈ Rny is the reference signal. To solve the optimization problem, and we used Levenberg-Marquardt Algorithm (Moré, 1978).
The result of the tracking control up to t = 400 is shown in figure 2. In the figure, the applied control input ut, the measured output yt (red line), reference signal rt (dotted black line), and the true water levels in the tanks xt,1, xt,2 (blue lines) are shown. In addition, the green lines in each plot show uft and y f t at t = 400, respectively. As can be seen from the figure, appropriate MPC is performed using the obtained state estimator and predictor. Also, the state estimate x̂t and the state of the system xt are compared in Appendix C.
Then, to support the interpretation discussed in section 4, figure 3 shows the visualization of the state estimator and predictor at the operating point at t = 50, where the change in water level is slow and nonlinearity is expected to be weak. The dashed lines in figures 3(a)-(c) and the bitmap on the right in figure 4(d) show the corresponding properties of the linearized system of (23) at the operating point. Since linear systems also have degrees of freedom in how they coordinate in state space, a transformation matrix Φ := P †x ∂P ∂x̂ ∈ R
nx×nx is used to align the coordinate system. From figures 3(a) and (b), the state estimator E obtained by the proposed method agrees well with the Kalman filter for the linearized system of the true system, and it can be confirmed that the proposed method yields almost optimal state estimators as expected for subspace methods in an ideal situation with sufficient data. It can also be seen that the predictor is consistent with the true system.
Moreover, figure 4 was obtained at the operating point at t = 380, where the nonlinearity is considered to be more strongly observed. Although the results in this case are similar to those in the linearized system to some extent, differences in the characteristics can be seen. For example, in
figure 4(d), the Toeplitz structure of linear systems is partially lost, indicating the existence of nonlinearity. In figures 4(a) and (b), it can be seen that an estimator that is less dependent on the output and more focused on the recent input than the Kalman filter for the linearized system is obtained.
To verify the validity of placing the bottleneck, we evaluated the performance with different widths of the bottleneck nx̂. All settings except nx̂ and the realized value of the stochastic noise wt, vt are set to be identical to the case of nx̂ = 2 and 60 trials were conducted with different realizations of the stochastic noise.
The results are summarized in figure 5. The loss for the validation dataset in figure 5(a) shows that
the performance does not degrade until the width of the bottleneck nx̂ goes below the actual state dimension nx = 2. Also, the performance of the MPC calculated by ∑400 t=hp+1
lt(ut, xt,2) does not show degradation until nx̂ < nx = 2 as shown in figure 5(b), and these results support the use of the bottleneck to provide interpretability.

heading: 6 CONCLUSION
text: In this paper, we proposed a new nonlinear subspace identification method that produces state estimator and predictor directly usable for nonlinear MPC. It is theoretically shown that by training a neural network with a bottleneck that performs regression on input-output sequences of a length that provides observability in both past and future directions, a state estimator and a predictor are formed across the bottleneck layer. The relationship between this method and the subspace identification method, as well as the interpretation of the obtained model using this relationship, is also explained, and its use is demonstrated through an example.
This research is still in its early stages, and there is much work to be done in the future. In order to demonstrate the applicability to higher-dimensional problems, an example of identifying dynamics from a simple video is included in Appendix D. However, verifying the performance in more practical tasks is a subject for future work. Also, a more precise condition for yielding a state estimator is desired. The current proof makes use of the assumption that zero-input responses from all states have been learned, and while the fact that zero-input responses constitute states is often used in this type of research (Verdult & Scherpen, 2004), more precise conditions would be useful for experiment design. It also appears that the technique of inverse optimal control (Ab Azar et al., 2020) can be applied to the interpretation of the obtained state estimator. Finally, since noise with high-order dynamics is ubiquitous (Bak et al., 1987) and the synthesis of state estimator for such a system is prone to failure, it is important to develop a mechanism to focus on useful states.

heading: D SUPPLEMENTARY EXAMPLE
text: To demonstrate the applicability to high dimensional systems, the proposed method is applied to the problem of estimating the dynamics of video discussed in (Beintema et al., 2021). Consider an environment where a nonlinear force acts a ball and we observe it with a camera of 25× 25 pixels. The dynamics of the ball is
p̈x = β
( 1
p2x − 1 (1− px)2
) − γṗx + kux (30a)
p̈y = β
( 1
p2y − 1 (1− py)2
) − γṗy + kuy (30b)
where β = 1/200, γ = 0.79, k = 1/4, and (px, py) is the position of the ball, which is not measurable. (ux, uy) is the external force whose sample value at time t are given as ut. This system
has a stable equilibrium point (px, py) = (0.5, 0.5) and we simulate this system from this point with zero velocity. For observation, the pixel intensity at a position X,Y of the video is calculated by nonlinear equation
pixel(X,Y ) = max ( 0, 1− (X − px) 2 + (Y − py)2
r2
) + v (31)
where r = 0.25 and v ∼ N (0, 0.2042) is noise term. The video image yt is sampled from this continuous image with 25× 25 grid from (0, 0) to (1, 1). This system has 2(= nu) inputs, 4(= nx) states, and 25× 25 = 625(= ny) outputs. We discretize this system by forward difference with the time step ∆t = 0.3 and simulate with random external force (ux, uy) uniformly distributed in [−1, 1] to generate video data. On training, we choose the hyperparameters of n̂x = 4,hp = hf = 20. E,P and f̂ have three, three, and two layers respectively, where each layer contains 64 ReLU neurons, the images are flattened into vectors with 625 elements, and all layers are fully connected. The training dataset consists of a sequence {(ut,yt)}100019t=−19 , and other settings are the same as in section 5. In the following, we describe the validation of the obtained model and its application to denoising. The data used in the following are obtained under the same settings as in training, except that the realizations of random input and noise are different.
D.1 MODELING OF DYNAMICS
To verify that the resulting model correctly models the dynamics of the video, we compare the video from the open-loop simulation of the obtained model with the video from the target system for the same input (external force).
Here, the output image from the model ŷt is computed recursively without using yt as
ŷt = [ Iny 0ny×ny(hf−1) ] P ( E(upt , ŷ p t ), [ ut
0nu(hf−1)
]) , (32)
where ŷpt is a sequence of past model outputs defined as in (2). Note that in this test, the ball is stopped at the center (0.5, 0.5) at t = 1, and ŷt for t = 1 − hp, . . . , 0 are set as the image at this state.
Figure 8 shows the actual images without noise (above) and images from simulation with the model (below). As shown in the figure, the images obtained from the model simulation are close to the
actual image without noise, indicating that the dynamics are correctly learned. A video version of figure 8 is provided in the supplementary material as actual-simulated.mp4.
D.2 NOISE REDUCTION
Next, we show that the model obtained by the proposed method can be used for denoising the video. Here the denoised image ŷt is computed from the images y p t and inputs u p t of the past hp frames
and the current input ut as
ŷt = [ Iny 0ny×ny(hf−1) ] P ( E(upt ,y p t ), [ ut
0nu(hf−1)
]) . (33)
Note that unlike the open-loop simulation (32), the state estimator E uses a noisy output ypt instead of recursively using ŷpt .
Figure 9 shows the result of denoising, The bottom images are the denoised images ŷt obtained from the trained model, and the middle images are the images obtained from the actual system yt. From the figure, it can be confirmed that the noise is removed by the model obtained by
the proposed method. A video version of figure 9 is provided in the supplementary material as actual-noisy-denoised.mp4.

Review summaries:
The authors are proposing to use DNN to approximate observers and predictors for nonlinear systems. No convergence analysis of the 'said' observers are given. The empirical results are rather weak.

This paper proposes a new system identification and model predictive control method based on neural networks for non-linear dynamical systems. The basic idea behind the proposed method comes from subspace identification for linear systems. In this sense, the proposed method can be regarded as a non-linear extension of subspace identification by deep neural networks.

The authors introduce a system identification methodology catered towards the model predictive approach. They claim that the main contribution is that the neural network, which is used, has a bottleneck layer and represents input-output dynamics. They also provide a connection to the subspace state-space system identification method for linear dynamical systems.