Control of Magnetic Manipulator Using Reinforcement Learning Based on Incrementally Adapted Local Linear Models

Research Article Control of Magnetic Manipulator Using Reinforcement Learning Based on Incrementally Adapted Local Linear Models Martin Brablc , 1 Jan ˇ Zegklitz , 2 Robert Grepl , 1 and Robert Babuˇ ska 2 , 3 1 Institute of Solid Mechanics, Mechatronics and Biomechanics, Faculty of Mechanical Engineering, Brno University of Technology, Brno 616 69, Czech Republic 2 Czech Institute of Informatics, Robotics and Cybernetics, Czech Technical University, Prague, Prague 16 636, Czech Republic 3 Cognitive Robotics, Delft University of Technology, Delft 2628 CD, Netherlands Correspondence should be addressed to Martin Brablc; [email protected] Received 11 December 2020; Revised 11 November 2021; Accepted 16 November 2021; Published 20 December 2021 Academic Editor: Aydin Azizi Copyright ©2021 Martin Brablc et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Reinforcement learning (RL) agents can learn to control a nonlinear system without using a model of the system. However, having a model brings benefits, mainly in terms of a reduced number of unsuccessful trials before achieving acceptable control performance. Several modelling approaches have been used in the RL domain, such as neural networks, local linear regression, or Gaussian processes. In this article, we focus on techniques that have not been used much so far: symbolic regression (SR), based on genetic programming and local modelling. Using measured data, symbolic regression yields a nonlinear, continuous-time analytic model. We benchmark two state-of-the-art methods, SNGP (single-node genetic programming) and MGGP (multigene genetic programming), against a standard incremental local regression method called RFWR (receptive field weighted regression). We have introduced modifications to the RFWR algorithm to better suit the low-dimensional continuous-time systems we are mostly dealing with. The benchmark is a nonlinear, dynamic magnetic manipulation system. The results show that using the RL framework and a suitable approximation method, it is possible to design a stable controller of such a complex system without the necessity of any haphazard learning. While all of the approximation methods were successful, MGGP achieved the best results at the cost of higher computational complexity. Index Terms–AI-based methods, local linear regression, nonlinear systems, magnetic manipulation, model learning for control, optimal control, reinforcement learning, symbolic regression. 1. Introduction A reinforcement learning (RL) agent interacts with the system to be controlled by measuring its states and applying actions according to a policy so that a given goal state is attained. The policy is iteratively adapted in such a way that the agent receives the highest possible cumulative reward, which is a scalar value accumulated over trajectories in the system’s state space. The reward associated with each transition in the state space is described by a predefined value function. Existing RL algorithms can be divided into critic-only, actor-only, and actor-critic variants. The critic-only variants optimize the value function (V-function) that is then used to derive the policy; the actor-only variants work directly on the policy optimization without any need for a value function; and actor-critic variants optimize both functions simultaneously. An example of the actor-only RL variant, often called Q-learning, can be found in [1] and that of the actor-critic variant in [2]. From a different point of view, RL algorithms can be also divided into model-based and model-free variants. Examples of both approaches can be found in [3, 4]. The modelbased variants include a model representation of the system to be controlled and can be pretrained in simulation (offline) and then updated when controlling the actual system (online). Model-free methods learn online exclusively through trial and error. Both variants have their specific advantages and disadvantages. We can often find remarks about the model-free approach requiring much more data, Hindawi Complexity Volume 2021, Article ID 6617309, 12 pages https://doi.org/10.1155/2021/6617309 especially in high-dimensional cases [1]. In this paper, we employ the model-based, critic-only variant without any online training so that we can compare different modelling approaches. We focus on two promising categories of approximation algorithms: genetic programming and local linear regression. Our aim is to contribute to the methodology of choosing the optimal out of dozens of existing modelling algorithms when presented with a specific RL task. This problem arises not only in connection with the RL framework (see [5]) but in modelling of a dynamical system in general [6, 7]. Genetic algorithms (GA) and their many variations are well established as a tool for modelling or parameter estimation of dynamical systems [8, 9]. However, genetic programming as a modelling approach used within RL is relatively new and promises good results with high-dimensional systems where other approaches fail. It creates a continuous-time, globally nonlinear model described by an analytical equation built of combinations of predefined functions [10]. As it is common with genetic optimization algorithms, these methods tend to be computationally demanding. On the other hand, local regression is a wellestablished modelling approach for model-based RF agents where the model is composed of local linear models, offering fast and computationally cheap approximation. There are several variants of local modelling methods; comprehensive examples of grid-based local linear model structure and data-based local linear regression (LLR) are described in [11, 12], respectively. Even though the use of local regression techniques within RL has been researched in the past, it was mainly based on simple, memory-based approximation methods such as the LLR, which is thoroughly described and examined in [13, 14], and more complex incremental methods such as the receptive field weighted regression (RFWR) [15, 16] or locally weighted projection regression (LWPR) [17] were omitted, with the exception of [18], where the RFWR algorithm was used as a critic approximator. The RFWR and LWPR methods provide significant benefits in lower memory use and higher stability by employing optimization-based (RFWR) or statistical (LWPR) methods to discover the optimal distribution of the local models’ areas of validity, that is, the receptive fields. It is important to benchmark the modelling methods because of the large number of existing approaches, which aim at similar tasks, while there are no simple guidelines on the method choice. Also, the presented algorithms are not yet well established within the RL domain. Finally, studying control algorithms for magnetic manipulation systems has importance on its own because of its application in many industrial fields (medical applications, magnetic levitation systems, etc.), thus leading to the two separate aims of this paper: exploring control algorithms suitable for control of precise magnetic manipulator systems and benchmarking different modelling approaches. When dealing with real magnetic manipulator systems, we also need to address practical issues that are often neglected in simulations, that is, nonlinearities such as actuator dead zones, saturations, Coulomb friction, signal delays, and so on. These present significant obstacles then implementing the control algorithm on a real system. In some cases, the dead zone and saturation problem can be addressed by nonlinear or adaptive control laws. For example, [19] shows an approach using fuzzy control with Gaussian membership functions, which is in practice similar to the RFWR method, and [20] describes a gain-scheduling adaptive approach to deal with internal system bounds. Using RL to find a control law for a nonlinear system also has the advantage that it can often deal with such disturbances on its own through the optimization process; for example, only a limited range of the actor outputs may be limited, which is the approach utilized in this paper. In this paper, we also present minor adjustments to the RFWR algorithm in Section 3, proposed to lower the computational complexity while preserving stability when working with low-dimensional problems. 2. Methods 2.1. Magnetic Manipulator. Genetic programming was already applied to nonlinear systems like an inverted pendulum or a collaborative robot [2, 10, 21]. To further investigate the approximation capabilities of these methods, we use a different system–a magnetic manipulator (Magman). This system consists of four coils that are independently operated by separate current controllers and a steel ball that can move freely over the coils; see Figure 1. To ensure that the ball moves only in the measured direction with limits on the edges, it is placed in a groove with 10 mm in size. In this case, we decided to limit the system to the first two coils only, as a system with four inputs is much more complex in terms of the RL computational complexity, while it does not enrich the system with different nonlinearities as it only spatially repeats the same kind of nonlinear behaviour. The steel ball can be positioned by properly controlling the current and thereby the magnetic force of the coils. The magnetic force a coil exerts on the ball is highly dependent on the distance of the ball from the coil’s centre, which introduces a significant nonmonotonic nonlinearity [22–24]. All experiments and simulation were scripted in MATLAB. The coil currents are controlled by stabilized current source modules, which communicate that MATLAB through a USB/RS232 transceiver using the virtual COM port (VCP) protocol on Windows OS. As the ball position is measured with a laser sensor with analog (voltage) output, the Humusoft MF634 IO card was used to measure the signal in real time from the MATLAB environment with a sampling period of 5 ms. Even though the Window OS is not an RTOS, with this sampling frequency, the period jitter is negligible (below 0.1%), and thus, the system can be considered real time. Table 1 lists the parameters of the magnetic manipulator we use in our experiments. With the task being a precise positioning of an object in a magnetic field, similar concepts can be found in many real-world applications, for example, maglev, microrobots, contactless stirring of chemicals, and so on. 2Complexity Approximate equations of motion inferred using the first principle method can be found in [24]. The system parameters were either measured directly or estimated using MATLAB parameter estimation toolbox based on measured data. Generally, the system can be described by a continuoustime, nonlinear state-space model as follows: x .�f(x, u),(1) where x� [x, _ x]⊤is the state vector composed of the position xand velocity _ xof the ball, forming the continuous system state space x∈X⊂R2;x .in _ X⊂R2is the state vector derivative; and u� [u1, u2]⊤is the input (action) vector composed of the coil currents. u1, u2form the system input space u∈U⊂R2. The nonlinear vector function f:X×U⟶_ Xthus describes the system dynamics. In this paper, by modelling the system, we mean approximating the underlying real function fusing various methods, which all build upon experimentally measured input-output data. Each data point is formed by corresponding assumed inputs and outputs of the function f –(x . k,xk,uk). In practice, these data points are corrupted by noise and other disturbances that are assumed to be with zero means. 2.2. SNGP. Single-node genetic programming (SNGP) is a graph-based genetic programming algorithm evolving a population organized as an ordered linear array of interlinked individuals, each representing a single program node [2, 10, 21]. Generally, symbolic regression algorithms try to find a model in the form of an analytical expression for a given data set by forming and evolving the expression out of elemental functions and operations. In our case, the algorithm is based on the assumption that the nonlinear function fin (1) can be efficiently approximated by the following equation: x .�f(x, u) � 􏽘 nf,ns i,j βifixj, uj 􏼐 􏼑,(2) where the nonlinear function fi, called the feature, is developed by means of genetic programming with nfbeing the maximum number of features, nsnumber of states, and the coefficients βiestimated by the least-squares method. The features are constructed from a list of elementary functions that are assumed to be able to produce the required fitting approximation of the presented data. The features can be combined by common operators or nested, but the maximal depth of the expression is limited to avoid overfitting. The symbolic model is evolved so that the mean-squared error over the training data is minimized. 2.3. MGGP. The second GP algorithm we used is called multigene genetic programming (MGGP). As opposed to SNGP, it combines the features defined also by 2 into treelike structured expressions called genes. The final expression is formed by a linear combination of these genes, which act as the individual features in equation (2). The parameters of this top-level linear combination are again estimated through least squares. Further details about the algorithm can be found in [25]. The actual MGGP implementation we used is extended with linear combinations of features [26] that enable the algorithm to find affine transformations of the feature space via a backpropagation-like technique, thus making it easier for the driving genetic programming algorithm to approximate the nonlinearities. 2.4. Receptive Field Weighted Regression. Receptive field weighted regression (RFWR) is an incremental approximation method that creates a set of local linear models and the corresponding Gaussian basis functions called the receptive fields and gradually updates them to fit the inputoutput data. The set of local linear models is updated with new data points (called the query points) using a weighted variant of the recursive least squares (RLS) method and the basis functions are updated through a gradient search with the help of heuristic decision rules. It can continually improve the set of models while still providing the best estimation of the approximated function at each query point based on the previously provided data. The original algorithm, first presented in [15, 16], which is the basis we build upon, can be best described by the following pseudocode: (1) For each new query point (x . k,xk,uk) (2) For each existing local model (3) Calculate model weight waccording to (4) (4) If w>activation limit wact (5) Update model parameters using RLS according to (6) and (7) (6) Update the corresponding receptive field using (12) and (14) (7) End (8) End Figure 1: A schematic drawing of the magnetic manipulator system. Table 1: Magnetic manipulator parameters. Parameter Value Unit Ball mass 53 g Ball diameter 20 mm Distance between edge coils and the ball position limit 20 mm Distance between coils 25 mm Maximal coil current 0.6 A Sampling period 0.005 s Complexity 3 (9) If no model was activated (10) Place a new model at the query point using (15) (11) Else if two or more models were activated with weight w>pruning limit wprun (12) Prune the model with the smaller receptive field (13) End (14) Calculate the model output as a weighted average of the activated local models (15) End Usually, the receptive field activation limit is set as wact �0.001. This parameter represents the weight limit for a local model to be updated according to the new data and to be included in the output estimation through a weighted average with another activated model. The pruning limit is usually set as wprun �0.7, which represents the highest acceptable overlap of neighbouring receptive fields. The RFWR variant described in this paper follows the main outline of the original algorithm [15] with several adjustments and improvements for the sake of stability and computational complexity for low-dimensional problems. This mainly concerns the rules for adding new local models, adjusting their receptive fields, and generalizing the algorithm in a way that the receptive fields are placed and optimized in a lower number of dimensions than the order of the models. This is especially useful in cases when the nonlinearities are significant mainly in one or two dimensions of the state space of the system. This algorithm, in its original implementation, is successfully being used to approximate inverse models of nonlinear systems to be used as a feedforward compensator [27, 28]. Figure 2 shows an example approximation of a complex univariate nonlinear function by the RFWR algorithm. Each of the local models is represented by a parameter vector b� [b1, b2,. . . , bn]⊤. With the input vector (a query point) Xq� [x1, x2,. . . , xn]⊤, the output yqis calculated by yq�X⊤ qb. (3) The weight wof a local model at a query point Xqis determined by its Gaussian receptive field as follows: w(X) � e−1/2 Xq−c  􏼁⊤C−1Xq−c  􏼁,(4) with c� [c1, c2,. . . , cn]the vector of model centre coordinates and C−1the distance inducing matrix of the basis function (receptive field). The overall output is then calculated as a weighted average of the outputs of the activated local models. The output estimate of the set of local models and their receptive fields is calculated by the following equation: 􏽢 y�1 􏽐n i�1wi􏽘 n i�1 wiyi.(5) We modified the original RFWR algorithm described in [15] to be used for low-dimensional problems. These modifications consist of the following: (1) Precise placement of new receptive fields that takes into account the location and dimensions of the existing surrounding receptive fields (see section 3.3) (2) Heuristic rules for stable updating of the receptive fields (see section 3.1) (3) Receptive fields can be distributed along a lower number of dimensions than the dimension of the data space (see section 3.2) 2.4.1. Updating parameters of Local Models. When a new query point is acquired, the activated local models are updated using the recursive least-squares algorithm (RLS) according to Pn+1�1 λPn−PnXT qXqPn λ/w+XqPnXT q ⎛ ⎝⎞ ⎠, bn+1�bn+wPn+1XT qe, e�yq−Xqb, (6) x y Nonlinear function approximation by RFWR Nonlinear function Local models x Receptive field weight Receptive fields Figure 2: Example of RFWR approximation of a complex nonlinear function and the distribution of receptive fields. 4Complexity where Pis the covariance matrix of the estimate, λis a forgetting parameter, and yqis the acquired output for the actual system state Xqcalled the query point. The covariance matrix Pneeds is usually initialized as a diagonal matrix. 2.4.2. Updating dimensions of Basis Functions. To avoid calculating the matrix inversion in (4) for every local model, an upper triangular matrix Mis used instead of C. Because of symmetry and positive definiteness, these matrices relate according to C−1�MTM. (7) To update the receptive field, we update Musing a gradient-descent optimization Mi+1�Mi−αzJ(M) zM,(8) of the cost function Jas follows: J�1 􏽐n i�1wi􏽘 n i�1 wiyq−yi 􏼐 􏼑2,(9) where wiis the activated receptive field weight, yiis the estimated output of the respective model at the query point (yq;Xq), and nis the number of local models. The parameter αis the gradient optimization step size. As the calculation of the cost function Jaccording to (9) is computationally very complex, we simplified the optimization algorithm through a set of heuristic decision rules and implemented the optimization as follows: Mi+1�Mi−αpzw(M) zM zw(M) zM� ze−1/2 Xq−c  􏼁⊤MTM Xq−c  􏼁 􏼒 􏼓 zM� � − Xq−c 􏼐 􏼑⊤M Xq−c 􏼐 􏼑e−1 2Xq−c 􏼐 􏼑⊤MTM Xq−c 􏼐 􏼑. (10) This implementation introduces a parameter p, which is an expression of a simple heuristic to decide whether the value of a basis function (weight) at the actual query point should be increased or decreased. This enables to stop updating the distance inducing matrix when a precision criterion is met and to limit the maximal number of local models to avoid overfitting. Parameter pcan be determined by various decision rules. A simple yet effective set, which was used in this research, can be created by using a long-term (cumulated over time) MSE of a particular model according to the data points, which can be described by p�−1,if MSE >MSElim 1,if MSE <MSElim 􏼨(11) 2.4.3. Adding New Local Models. During the optimization process, it is possible that no model exceeds the activation limit wact. In such a case, a new local model with a receptive field is added to the approximation set. The centre of the receptive field is automatically placed at the actual query point, and the model parameters are initialized to fit the measured output of the approximated system. What needs to be determined is the area in the state space that should be covered by the newly created receptive field. The original algorithm uses a default diagonal distance inducing matrix for every local model. However, an optimal distance inducing matrix can be determined. Intuitively, the new receptive field should cover the gap between the already existing models. The distance inducing matrix should be initialized as a diagonal matrix with parameters that ensure that the new receptive field does not overlap with any existing one more than a preset limit. In our case, the limit was set to 0.5wprun. Since this would be a complex optimization task not suitable for realtime calculation, we simplified the criterion so that the maximal overlapping weight of two models is analyzed only over the line segment connecting their centres. In that case, the distance parameter for initializing the distance inducing matrix can be determined by the following equation.where vi�cn−ciis a vector between the new centre cnand the centre of a neighbouring receptive field ci. A two-dimensional example is shown in Figure 3. This method yields a better estimate of the distance inducing matrix of the new receptive field than the fixed initial dimension matrix in the original algorithm as it requires fewer iterations to stabilize and to cover the gap between neighbouring receptive fields. di� − 2 log wprun/2 􏼐 􏼑 vi 􏼌􏼌􏼌􏼌􏼌􏼌􏼌􏼌−���������������������� −2 log wprun/2 􏼐 􏼑/viMT iMivi 􏽱,(12) The distance parameter dihas to be calculated for every existing local model, and the minimal distance dmin is used to initialize the distance inducing matrix according to M0 n�I���� dmin 􏽱,(13) where Iis a unity matrix of the corresponding order. In the specific case of the magnetic manipulator, the inputs of the local models would correspond to (x, u)and the output to x .. 2.5.ReinforcementLearning. Consider the following discrete deterministic state-space model of a system to be controlled: xk+1�f xk, uk  􏼁,(14) where k∈Zdenotes discrete time instants, xk, xk+1∈X⊂Rnis the state vector, and uk∈U⊂Rmis the input vector. An RL agent learns to control the system so that it achieves the maximal cumulated reward on a trajectory from the initial state to the desired state [10]. At each state transition, as described by (14), the agent receives a scalar reward according to Complexity 5 r�ρxk, uk, xk−1  􏼁.(15) The reward function ρis usually based on the distance of the current state to the goal state. The optimal control law, called the policy, π:X⟶Uis determined as follows such that it maximizes the cumulative reward, called the return: Rπ�E􏽘 ∞ k�0 ckρxk,πxk  􏼁, xk+1  􏼁 ⎧ ⎨ ⎩⎫ ⎬ ⎭,(16) where c∈(0,1)is the discount factor and the initial state x0 is selected from the state space domain X. The return for any permissible initial state xis captured by the value function V:X⟶Rdefined as follows: V(x) � E􏽘 ∞ k�0 ckρxk,πxk  􏼁, xk+1  􏼁 ⎧ ⎨ ⎩⎫ ⎬ ⎭, x0�x. (17) An approximation of the optimal V-function 􏽢 V(x)can be found by solving the Bellman equation as follows: 􏽢 V(x) � maxu∈U[ρ(x, π(x), f(x, u)) + c􏽢 V(f(x, u))]. (18) The optimal action can be found as the action that steers the system to a state with maximal value [21]. This corresponds to maximization of the right-hand side of (18): u�arg max u′∈U ρx, u′, f x, u′  􏼁 􏼁+cV f x, u′  􏼁 􏼁􏼂 􏼃.(19) 3. Experimental Results We prepared training and validation I/O data sets measured on the magnetic manipulator with random input signals as a list of data points in the form (x . k, xk, uk). The random input signals (coil currents) were generated in the way that only one coil was active at a time that eliminated possible electromagnetic interactions between them (the coil current was controlled by an HW-based current feedback controller module rendering the transient times negligible). Figure 4 shows an example of a training data set. Even though the ball’s position measurement is very precise, it still contains significant noise. For that reason, the time-domain derivatives of the position (velocity and acceleration) needed for the dynamic model approximation were determined using the Savitzky–Golay filter, which is an FIR filter based on least-squares polynomial approximation able to perform numerical differentiation while filtering the noise simultaneously [29, 30]. Especially for the RFWR implementation, it is important to note that the system’s nonlinearity is mainly significant along the position of the ball and the system can be seen as linear in parameters along the other dimensions (acceleration and velocity). In this case, the general model (1) can be rearranged as follows: € x�f(x, u) + b0sign(_ x) + b1_ x+b2_ x2,(20) where the function f(x, u)represents the significant nonlinearity suitable for local approximation, the term b0sign(_ x)represents a simple model of dry friction, the term b1_ xrepresents the viscous friction, and the last term b2_ x2 models nonlinear damping caused by electromagnetic induction influencing the steel ball while moving rapidly through a magnetic field. Despite being nonlinear, all of the terms are linear in their parameters and can be modelled globally, which means that the local models share parameters b0through b2. The term b0sign(_ x)in the (20) is quite important in practical situations where Coulomb friction is not negligible. The sign function is often being used to approximate the effects of Coulomb friction wherever there is no significant stiction (difference between static and dynamic friction effects). There are better approximations for simulation purposes, for example, the sigmoid function; however, most of them are not linear in parameters and thus not applicable for RLS parameter estimation. The same data set was presented to all of the approximation methods (RFWR, SNGP, and MGGP). Due to the stochastic nature of the two algorithms based on genetic programming, the same process was repeated with different pseudorandom seeds. Overall, 30 runs for SNGP and MGGP and 1 run for RFWR were made. Table 2 shows the summary of the MSE results. Since the MSE of the models with respect to the training set is not sufficient to decide which models are better, two Receptive fields - 2D Receptive fields - 1D section Old receptive field New receptive field X2 X1 ci ci wmax w cn X cn SECTION PLANE Figure 3: Demonstration of adding a new reference field to the set of RFWR local models and matching its dimensions to fit the surrounding reference fields. 6Complexity separate data sets were measured: one was used for the training of the models and one for validation. However, since the magnetic manipulator is not open-loop stable, the common open-loop validation is not suitable, as every model diverges quickly even though the parameters may be close to ideal due to errors introduced by numerical integration. Therefore, the models were validated in several steps-ahead prediction mode. As we also suspected that only onesample-ahead validation could be influenced by the remaining noise in the measured signals, we validated the model for 1, 3, 5, 10, 50, 100, and 250 samples ahead. We used the five-step-ahead prediction (5-SAP) as a baseline for selecting the best models for further experiments. The reason to choose five samples is based on an experimentally 0 0.5 t (s) 1 1.5 2 0 0.5 t (s) 1 1.5 2 0 20 40 60 80 x (mm) U (A) Measured 0 0.2 0.4 0.6 Coil 1 Coil 2 Figure 4: Random input signal (generalized binary noise with random coil switching) with only one coil active at a time for training and validation sets generation. Table 2: Statistical results of the best models. Model Training MSE 5-SAP validation MSE Control MSE SNGP 3 6.53 ·10−25.71 ·10−12 4.02 ·10−5 SNGP 7 6.78 ·10−25.77 ·10−12 2.46 ·10−4 SNGP 10 5.89 ·10−25.36 ·10−12 4.23 ·10−5 SNGP 13 7.28 ·10−26.37 ·10−12 4.28 ·10−5 SNGP 17 7.48 ·10−25.84 ·10−12 1.28 ·10−4 SNGP 18 7.23 ·10−25.79 ·10−12 2.06 ·10−4 SNGP 21 7.65 ·10−26.41 ·10−12 3.40 ·10−3 SNGP 27 5.32 ·10−26.57 ·10−12 1.16 ·10−4 SNGP 29 5.85 ·10−26.19 ·10−12 9.60 ·10−5 SNGP 30 6.93 ·10−25.99 ·10−12 3.75 ·10−5 MGGP 5 5.94 ·10−25.49 ·10−12 6.34 ·10−5 MGGP 8 4.87 ·10−25.13 ·10−12 4.14 ·10−5 MGGP 9 5.69 ·10−26.18 ·10−12 3.80 ·10−5 MGGP 13 6.59 ·10−25.67 ·10−12 4.13 ·10−5 MGGP 16 5.96 ·10−25.29 ·10−12 6.59 ·10−5 MGGP 17 5.99 ·10−25.08 ·10−12 5.46 ·10−5 MGGP 18 7.35 ·10−26.28 ·10−12 1.50 ·10−4 MGGP 21 5.60 ·10−25.92 ·10−12 6.74 ·10−5 MGGP 22 5.55 ·10−25.74 ·10−12 8.84 ·10−5 MGGP 23 6.00 ·10−25.82 ·10−12 3.56 ·10−5 RFWR 6.25 ·10−21.24 ·10−11 6.57 ·10−5 Complexity 7 validated assumption that shorter intervals do not show the model’s imprecision and longer intervals cause even very precise models to diverge randomly. The n-step-ahead prediction validation is based on a moving frame of nconsecutive data points, where the first data point is applied as the initial condition for numerical integration (using the ode45 solver) of the dynamical model being tested. When the simulation reaches the n-th step, an MSE residual is calculated between the corresponding data point and the model prediction. The resulting model validation metric is then calculated as the sum of residuals over each of the prediction frames. Figure 5 shows examples of the models constructed by each algorithm. As the transition model of the system is four-dimensional, for visualization purposes, the figures show a two-dimensional situation for the input vector set to u� [0,0.6]T. This corresponds to the situation when the first coil is turned off and the current through the second coil is 0.2 0 -4 SNGP model -3 -2 0 -1 0.02 x (m) 0 -0.2 1 0.04 2 0.06 3 0.08 00.02 x (m) x (m) 0.04 0.06 0.08 0.2 0 MGGP model -0.2 0.2 0 RFWR model 00.02 -0.2 0.04 0.06 0.08 x ·· (m·s-2) -4 -3 -2 -1 0 1 2 3 x ·· (m·s-2) -4 -3 -2 -1 0 1 2 3 x ·· (m·s-2) x · (m·s-1) x · (m·s-1) x · (m·s-1) Figure 5: Comparison of approximated dynamic models for the SNGP, MGGP, and RFWR method, best-fitting models from each category. The 3D visualization is shown for constant input u� [0,0.6]T. x (m) 0.1 -25 -20 -15 V (–) V (–) V (–) SNGP model 0 -10 0.02 -5 -0.1 0 0.04 x (m) x (m) 0.02 0.04 0.1 -25 -20 -15 MGGP model 0 -10 -5 -0.1 0 0.1 -25 -20 -15 RFWR model 0 -10 0.02 -5 -0.1 0 0.04 x · (m·s -1 ) x · (m·s -1 ) x · (m·s -1 ) Figure 6: Comparison of V-functions based on the best fitting models from each of the approximation methods: SNGP, MGGP, and RFWR. 8Complexity 0.02 0.04 -0.1 0 0.1 SNGP - Coil 1 SNGP - Coil 2 MGGP - Coil 1 MGGP - Coil 2 RFWR - Coil 1 RFWR - Coil 2 x · (m·s-1) -0.1 0 0.1 x · (m·s-1) -0.1 0 0.1 x · (m·s-1) -0.1 0 0.1 x · (m·s-1) -0.1 0 0.1 x · (m·s-1) -0.1 0 0.1 x · (m·s-1) x (m) 0.02 0.04 x (m) 0.02 0.04 x (m) 0.02 0.04 x (m) 0.02 0.04 x (m) 0.02 0.04 x (m) Figure 7: Comparison of policies based on the best fitting models of the magnetic manipulator from each of the approximation methods. 05 t (s) x (mm) U (A) 10 15 05 t (s) 10 15 0 10 20 30 40 0 0.1 0.2 0.3 0.4 Coil 1 Coil 2 Measured Desired Figure 8: Example of the control performance for a controller based on the SNGP 10 model. Complexity 9