# Feedforward Nonlinear Control Using Neural Gas Network

Authors: Iván Machón-González; Hilario López-García
Journal: Complexity
 (2017)

---

## Abstract

Nonlinear systems control is a main issue in control theory. Many developed applications suffer from a mathematical foundation not as general as the theory of linear systems. This paper proposes a control strategy of nonlinear systems with unknown dynamics by means of a set of local linear models obtained by a supervised neural gas network. The proposed approach takes advantage of the neural gas feature by which the algorithm yields a very robust clustering procedure. The direct model of the plant constitutes a piece-wise linear approximation of the nonlinear system and each neuron represents a local linear model for which a linear controller is designed. The neural gas model works as an observer and a controller at the same time. A state feedback control is implemented by estimation of the state variables based on the local transfer function that was provided by the local linear model. The gradient vectors obtained by the supervised neural gas algorithm provide a robust procedure for feedforward nonlinear control, that is, supposing the inexistence of disturbances.


---

## Body

## 1. Introduction

Although some physical systems can be approximated as a linear model, almost all real plants actually have a nonlinear functioning. A wide understanding of the behavior of nonlinear processes is available but it is sometimes difficult to choose the appropriate control method. Lyapunov theory is a classic method for nonlinear system control. If and only if there is a positive definite continuous function whose derivative is negative under the suitable conditions of the control design, then the control asymptotic stability is guaranteed. However, this method is unfortunate because obtaining the Lyapunov function is difficult. This problem is even worse when dealing with unknown plants that are not defined mathematically. Therefore, it is usually not easy to guarantee the stability of a complex nonlinear system [1]. However, if the local linear system corresponding to an equilibrium point is controllable, then sufficient conditions can be stated for local stability [2].Hartman-Grobman theorem states that the behavior of a nonlinear system in the neighborhood of an equilibrium point can be approximated by its linearized model. The systems theory is based on many mathematical procedures about stability, controllability, and observability regarding linear systems. The stability and, to a great extent, the dynamic response of a linear system can be described in terms of eigenvalues of the system matrix in state space design or poles of the transfer function. No such method exists for nonlinear systems. For this reason, industrial control processes are still usually designed using this linear control theory. After linearization, the typical approach is to design a linear controller such as PID with fixed parameters.The classical approach to get local linear models can be achieved with RLS (Recursive Least Squares) method. However, sometimes this method throws up unfavorable results due to the intrinsic nonlinearities of the process to be controlled. The problem is to establish the different operating points for a nonlinear system. At this point, the proposed algorithm can establish each operating point as a cluster centre of the neural gas network. It is for such reason that artificial intelligence techniques improve the control performance.Research into identification and control of nonlinear systems by means of neural networks (NN) began over two decades ago [3]. One of the major advantages of control by NN is that precise knowledge of the plant such as a mathematical model is not needed. Initially, control applications using NN were based on a trial-and-error approach. Research efforts have improved the control algorithms and several journals have published special issues with a strong mathematical foundation [4]. Many applications are based on a combination between feedforward and recurrent NN. Recurrency, also known as dynamic backpropagation, is necessary due to the dependency of the output on the previous values of the same output, which are also functions of the weights [5]. Zhang and Wang [6] proposed a pole assignment control using recurrent NN.The typical design procedure is to carry out the system identification in order to model the plant and, secondly, to obtain the controller. Traditional methods rely heavily on models extracted from physical principles, whereas approaches based on NN theory usually create black-box models as function approximators using data obtained from the plant. Knowledge about the mathematical model of the plant or any other physical principle is not necessary.Neural gas (NG) is an unsupervised prototype-based method [7] in which the prototype vectors are the weights and carry out a partition of the training data space. It considers the cooperation-competition computation, allowing the algorithm to be prevented from local minima problems. In addition, the batch NG allows fast training so that the convergence is achieved in a small number of epochs [8]. Supervised versions of NG have also been developed, specially for classification [9, 10]. The algorithm has a great robustness for clustering tasks but has also been proven to be robust to obtain direct models of plants [11].After years of works in identification and control of dynamical systems by means of NN, there is agreement among researchers that linear identifiers and controllers should be used as first attempt, as stated in Chen and Narendra [2]. If a set of local linear models corresponding to several equilibrium points can approximate with certain accuracy a nonlinear system, then linear controllers can be designed for each model and the global control is related to control by switched linear models.This divide-and-conquer approach is applied in this work. The resulting model is a set of local linear maps. Each neuron of the NG model corresponds to one local model. These local models are obtained after NG training. In this way, a direct model of the plant is obtained. After obtaining this NG model, the design of the local linear controller is simpler than that of the global nonlinear controller. Local linear mapping using another prototype-based algorithm such as SOM was successfully tested at the NASA facilities [12].This paper aims to apply the robustness modeling capability of NG to control a nonlinear plant such as a typical robot manipulator.The paper contains the learning rules of the considered NG algorithm in Section2, the model of the plant and the control strategy are explained in Sections 3 and 4, respectively, and the proposed technique is tested in Section 5.


## 2. Neural Gas Approach

The unsupervised version of the NG algorithm is based on energy cost function (1) according to the Euclidean metric. The notation used for the squared Euclidean distance is given in (2). Moreover,(1)ENG=∑i=1m∑j=1Nhσvj,wi·dvj,wi,(2)dvj,wi=vj-wi=vj-wi2.A neighborhood function (3) is needed to implement the algorithm. The rank function k(v,wi)∈0,…,m-1 represents the rank distance between prototype wi and data vector v. The minimum distance takes the value 0 and the rank for the maximum distance is equal to m-1, where m is the number of neurons or prototypes and σ(t) is the neighborhood radius:(3)hσv,wi=exp-kv,wiσt.The neighborhood radiusσ(t) is usually chosen to decrease exponentially according to (4). The decrease goes from an initial positive value, σt0, to a smaller final positive value, σtmax:(4)σt=σt0·σtmaxσt0t/tmax,where t is the epoch step, tmax is the maximum number of epochs, and σt0 was chosen as half the number of map units (σt0=m/2), as in Arnonkijpanich et al. [13]. In addition, σtmax=0.0001 in order to minimize the quantization error at the end of the training.The learning rule of the batch version is obtained in Cottrell et al. [8]. The batch algorithm can be obtained by means of Newton’s method using the Jacobian and Hessian matrices, J and H, respectively, of the cost function ENG. The adaptation of the prototype wi is formulated accordingly based on this method(5)Δwi=-Jwi·H-1wi.Kernel functionhNG can be considered locally constant [8]. In this way, the Jacobian and Hessian matrices are(6)Jwi=-2·∑j=1Nhσvj,wi·vj-wi,Hwi=2·∑j=1Nhσvj,wi.Substituting (6) into (5), the increment can be obtained(7)Δwi=∑j=1Nhσvj,wi·vj-wi∑j=1Nhσvj,wi.Finally, the updating rule for each prototype vector appears in(8)wi=∑j=1Nhσvj,wi·vj∑j=1Nhσvj,wi.


### 2.1. Supervised Learning

Supervised learning with NG is possible by means of local linear mapping over each Voronoi region defined by prototype vectorwi. A constant yi and a vector ∇fi with the same dimension as wi are assigned to each neuron i. The goal is to approximate the function y=f(v) from RD to R, where D is the number of training variables, that is, the dimension of data vector v. The training thus becomes supervised and the dataset contains input-output pairs of data vector v and variable y as the objective function. The estimation is carried out by (9)y^v=yi∗+∇fi∗·v-wi∗,where y^(v) is the estimated output value, yi is the reference value learned for wi, ∇fi is the gradient of the approximated function obtained in the ith Voronoi region defined by wi, and i∗ is the neuron i with its closest wi to data vector v, that is, the best matching unit (BMU). The asterisk super index denotes the winning neuron for input data vector v.The probability distribution of the input data is represented by prototype vectorsw which are previously updated according to the typical rule of the unsupervised version of the algorithm [14] using (8). Each prototype vector wi can be considered as the centroid of the ith Voronoi region. After unsupervised training, m regions are well defined by these vectors. At this point, local models will be created in each of these regions so that m local models will represent the whole data distribution.The energy cost function of the supervised version of the algorithm is based on the mean square error of the output variable estimation averaged over each Voronoi region [15] according to (10). Prototypes wi are already obtained in (8), whereas the adaptation rules of yi and ∇fi are calculated considering Newton’s method for energy cost (10). The learning rules for yi and ∇fi are shown in (11) and (12), respectively:(10)ENGsup=∑i=1m∑j=1Nhσvj,wi·fvj-f^vj2=∑i=1m∑j=1Nhσvj,wi·yj-yi-∇fi·v-wi2,(11)yi=∑j=1Nhσvj,wi·fvj∑j=1Nhσvj,wi,(12)Δ∇fi=∑j=1Nhσvj,wi·vj-wi·fvj-yi-∇fi·vj-wi∑j=1Nhσvj,wi·vj-wi·vj-wi.


## 2.1. Supervised Learning

Supervised learning with NG is possible by means of local linear mapping over each Voronoi region defined by prototype vectorwi. A constant yi and a vector ∇fi with the same dimension as wi are assigned to each neuron i. The goal is to approximate the function y=f(v) from RD to R, where D is the number of training variables, that is, the dimension of data vector v. The training thus becomes supervised and the dataset contains input-output pairs of data vector v and variable y as the objective function. The estimation is carried out by (9)y^v=yi∗+∇fi∗·v-wi∗,where y^(v) is the estimated output value, yi is the reference value learned for wi, ∇fi is the gradient of the approximated function obtained in the ith Voronoi region defined by wi, and i∗ is the neuron i with its closest wi to data vector v, that is, the best matching unit (BMU). The asterisk super index denotes the winning neuron for input data vector v.The probability distribution of the input data is represented by prototype vectorsw which are previously updated according to the typical rule of the unsupervised version of the algorithm [14] using (8). Each prototype vector wi can be considered as the centroid of the ith Voronoi region. After unsupervised training, m regions are well defined by these vectors. At this point, local models will be created in each of these regions so that m local models will represent the whole data distribution.The energy cost function of the supervised version of the algorithm is based on the mean square error of the output variable estimation averaged over each Voronoi region [15] according to (10). Prototypes wi are already obtained in (8), whereas the adaptation rules of yi and ∇fi are calculated considering Newton’s method for energy cost (10). The learning rules for yi and ∇fi are shown in (11) and (12), respectively:(10)ENGsup=∑i=1m∑j=1Nhσvj,wi·fvj-f^vj2=∑i=1m∑j=1Nhσvj,wi·yj-yi-∇fi·v-wi2,(11)yi=∑j=1Nhσvj,wi·fvj∑j=1Nhσvj,wi,(12)Δ∇fi=∑j=1Nhσvj,wi·vj-wi·fvj-yi-∇fi·vj-wi∑j=1Nhσvj,wi·vj-wi·vj-wi.


## 3. Plant Model

After NG training, the plant is modeled as a set of linear systems whose outputY depends on the previous values of both output Y and input U. The Nonlinear Autoregressive-Moving Average (NARMA) model has been proven for nonlinear identification [16, 17] and can be expressed as(13)Yk+d=hYk,Yk-1,Yk-2,…,Yk-n+1,Uk,Uk-1,Uk-2,…,Uk-n+1,where Yk is the system output at the sampling instant k, Uk is the system input at instant k, and d is the system delay.Considering zero delay system and substituting (13) for (9) remains(14)Yk=yi∗+∇fi∗,1Yk-1-wi∗,1+∇fi∗,2Yk-2-wi∗,2+⋯+∇fi∗,nYk-n-wi∗,n+∇fi∗,n+1Uk-wi∗,n+1+∇fi∗,n+2Uk-1-wi∗,n+2+⋯+∇fi∗,2n+1Uk-n-wi∗,2n+1.Hereafter, the gradients will be denoted as coefficientsai and bi.∇ f i ∗ , 1 = a 1, ∇fi∗,2=a2,…,∇fi∗,n=an, ∇fi∗,n+1=b0, ∇fi∗,n+2=b1,…,∇fi∗,2n+1=bn.And the following terms will be gathered to form variableη:(15)ηk=yi∗-∑j=12n+1∇fi∗,j·wi∗,j.Denoting the polynomials with backward shift operatorz-1 by Az-1=1-a1z-1-a2z-2-⋯-anz-n and B(z-1)=b0+b1z-1+b2z-2+⋯+bnz-n(16)Yz=Bz-1Az-1Uz+1Az-1ηzwhich reminds one of an ARMAX model, η is not only a zero mean independent identically distributed white noise process but also a known disturbance calculated according to (15), and it depends on the input and output values since it is obtained by BMU i∗. The internal noise of the system can be included in η.Using thez-transform, Y(z) is the system output and U(z) is the system input where the controller must be connected.


## 4. Local Linear Control by State Feedback

If the system is linear (locally), then the superposition theorem can be applied remaining the linear transfer function between the system output and the control input as follows:(17)Yz=b0+b1z-1+b2z-2+⋯+bnz-n1-a1z-1-a2z-2-⋯-anz-nUz.Define(18)Qz=Uz1-a1z-1-a2z-2-⋯-anz-nand choose the following relationship between the state variables: x1(z)=z-nQ(z), x2z=z-n+1Qz,…,xn(z)=z-1Q(z).Transfer function (17) can be expressed in control canonical form for linear state space design as(19)x1k+1x2k+1⋮xn-1k+1xnk+1=010⋯0001⋯0⋮⋮⋮⋱⋮000⋯1anan-1an-2⋯a1x1kx2k⋮xn-1kxnk+00⋮01uk,(20)Yk=bn+anb0bn-1+an-1b0⋯b1+a1b0x1kx2k⋮xnk+b0uk.Assuming that the system is controllable, the purpose of the control by state feedback via pole placement is to assign a set of pole locations for the closed-loop system that will correspond to satisfactory dynamic response in terms of rise time, settling time, and overshoot of the transient response. The control law is a linear combination of the state variablesxi which are estimated in (19) by way of local transfer function (17).The characteristic polynomial of the closed-loop system depending on system matrixF, input matrix G, and gain vector K is(21)det⁡zI-F-GK=0whereas the characteristic polynomial of the desired pole locations is(22)pz=∏j=1nz-λj=zn+p1zn-1+p2zn-2+p3zn-3+⋯+pn.For annth-order system, the gain vector K=K1K2⋯Kn for state feedback is obtained by matching coefficients in (21) and (22) forcing the closed-loop poles to be placed at the desired locations:(23)K=p1+a1p2+a2⋯pn+an.It is possible that there are enough degrees of freedom to choose arbitrarily any desired root location by selecting the proper valuesKi. This is an inexact procedure that may require some iteration by the designer. The solution of the local linear model lies in finding the matrix or the regulator coefficients that implement the state feedback control. The stability condition for linear discrete-time systems is that all the eigenvalues must be inside the unit circle. Obviously, this criterion is not valid for nonlinear systems but there is a region inside the stable linear area where the asymptotic stability of the switched linear systems is achieved [18]. Thus, not only a desired dynamic response can be designed, but also stability criteria will be accomplished. In this work, this stability region was found by means of trial-and-error with different eigenvalues.The proposed control strategy scheme is shown in Figure1. Gain vector K is calculated to fulfill the dynamics according to (22) depending on the local linear model defined by the current winning neuron i∗ or BMU. State variables xi are also obtained by the local linear model of the NG in (19). Tracking of the setpoint reference is possible using the inverse static gain of the feedback loop. In addition, since disturbance η is known (it is included by the model), it can be compensated as -η/Ai∗(z-1). The transfer function of the prefilter has been chosen as (1-λprefilter)n/(z-λprefilter)n and determines the switching rate of the local linear models. Although the pole assignment method does not affect the zeros of the plant, the prefilter can be optionally designed in order to cancel dominant zeros located inside the unit circle.Figure 1
Control strategy by state feedback and local linear models.


## 5. Experimental Testing

The aim is to control the typical robot arm problem depicted in Figure2. Hagan et al. [19] focused on this plant to be controlled by dynamic propagation algorithm using a Model Reference Control architecture. Obviously the proposal control is not based on the mathematical model of the plant, but this is a well-known second-order nonlinear differential equation:(24)ut=Jθt¨+Bθt˙+MgLsin⁡θt,where M = 1 kg, J = 1 kgm2, L = 1 m, B = 1 kgm2/s, and g = 10 m/s2. The viscous friction coefficient is important regarding stability and dynamic response. If B=0, then the system is unstable in open loop and the necessary training data can not be obtained by open loop simulation of the plant. In the present work, the system is simulated in open loop to acquire the training data considering B=1 in order to propose a plant with more oscillatory response in comparison to Hagan et al. [19] where B=2. Plant input u is a uniformly distributed random signal with −4 and 8.1 as minimum and maximum values of amplitude, respectively. The pulse width must be carefully selected in order to model the transient and steady states correctly. Thus, the pulse width is equal to 14 seconds. Since NG is a vector quantization algorithm where the neurons are updated according to the probability distribution function of the training data, it is important to obtain a uniform distribution of output value θ in the training data. After plant simulation, it was observed that the present system output θk did not depend on the present system input uk in the training dataset. If n=2 is considered, then the plant is not correctly modeled, only the steady state is approximated but not the transient. If n=3, then the plant model is quite accurate but the control is not well performed due to the system nonlinearities. The optimum value is n=4. Batch NG training was carried out considering 50 epochs and 49 neurons to obtain the direct model of the plant. Fewer neurons cause a similar effect to that mentioned above for n=2 and using more neurons does not improve the control. Figures 3 and 4 show the results for testing data after training for n=4. Obviously, the sampling time must meet the criteria of Nyquist-Shannon theorem. The considered sampling time was 0.1 seconds.Figure 2
Plant used to test the proposed control.Figure 3
Results for testing data usingn=4.Figure 4
Results for testing data usingn=4.Eigenvaluesλj in (22) provide asymptotic stability and suitable dynamic response. We consider that all the eigenvalues are equal; that is, p(z)=(z-λ)4. In order to tune λ, the control system was tested for different eigenvalues and prefilter poles. Steps of amplitude 0.9 at the reference setpoint were used since θ values close to 0.9 represent the worst system working zone with high nonlinearity. Obviously, θ values over 0.9 are even worse, but we refer to the training data range. Eigenvalue λ was incremented from 0.2 to 0.9 in series of 0.1 amplitude steps. Figures 5 and 6 show the results. If the prefilter has a wide bandwidth (λprefilter is low), then lower values of λ produce instability. In Figure 5, lower λ values yield a considerable effect of the modeled disturbance η and the overshoot is high, whereas there are some rebounds between two BMU linear models for higher λ values. The lower λ, the lower the rise time (the wider the bandwidth). The optimum λ range is [0.5–0.6] for λprefilter=0.4. Thus, when using switched linear models there is a stability region of λ within the global stability area of the linear systems theory (unit circle) depending on the switching rate of the NG local linear models [18]. Here the switching rate is determined by the prefilter transfer function. If λprefilter is increased then good results are obtained for λ=0.2, and see Figure 6.Figure 5
Results forλprefilter=0.4.Figure 6
Results forλprefilter=0.9.Once the parameters had been determined, the NG approach was tested to track a constant reference for control of position in Figure7 and a variable reference such as a sinusoidal signal to check control of velocity in Figure 8. The linear estimation output θ^L is not the NG estimation θ^ in (9) but it is calculated by means of (20) and adding the value of known disturbance η.Figure 7
Results forλprefilter=0.7 and λ=0.2.Figure 8
Results forλprefilter=0.7 and λ=0.2.The worst control is when the reference value is close to 0.9. To illustrate this problem, PI control with fixed parameters was compared to NG approach. Obviously, PID control would achieve a fast and stable response due to derivative action but at the cost of an unrealizable control action and the system control would be affected by noisy signals. The linearized model using thes-transform in the neighborhood of θ0, denoted as Gθ0, is(25)ΘsUs=Gθ0=1Js2+Bs+MgLcos⁡θ0.A suitable PI design regarding settling time and smoothness response isPI(s)=kp·s+(ki/kp)/s=0.3·s+10/s considering θ0=0.5. Figure 9 shows that the system controlled by this PI with fixed parameters becomes more oscillatory when considering θ0=0.9 because the two complex poles of the plant change their location involving the change of the closed-loop poles. The two complex closed-loop poles are more dominant and, therefore, the system increases oscillation. In this design the trade-off is between the settling time and the oscillatory component of the response. The tuning is to change kp keeping constant the location of the zero in −10; that is, ki/kp=10. Figure 10 shows the influence of the adjustment of kp in the control design and the NG approach is displayed to be compared to that one.Figure 9
PI control at two different linearization points.Figure 10
NG in comparison to PI.The control strategy described above is valid for feedforward control, that is, supposing the inexistence of disturbances. In these conditions, the algorithm promises very good performance. However, disturbance rejection can be achieved adding an extra state variable so thatx˙(t)=e(t), where e is the tracking error, and following the steps described from (21) to (23), Figure 11 shows the rejection of a constant output disturbance of amplitude 0.05.Figure 11
Disturbance rejection.


## 6. Conclusions

In this paper a supervised version of neural gas (NG) algorithm is proposed to control nonlinear systems whose dynamic mathematical models are unknown. The identification of the plant is achieved with the NG model. In comparison to other types of neural networks, the formation of the NG model is a robust procedure since there are neither problems of local minima nor overfitting. The training data must be carefully selected in order to model the transient and steady states correctly. The NG algorithm tends to model the steady states quite well. Obviously, the transient must be correctly modeled in order to control the plant. In this way, the number of delayed samplesn and the number of neurons m are key parameters. There must be a sufficient number of neurons but the control is not improved if it is too large.The trained NG network produces a set of piece-wise local linear models. Each of these is represented by a neuron. The global controller is a set of linear controllers which are obtained by state feedback via pole assignment. This control does not affect zeros but if these are inside the unit circle, then they can be cancelled by the poles of the prefilter.Eigenvaluesλ inside the unit circle do not guarantee the asymptotic stability because the plant to be controlled is nonlinear. Therefore, the stability corresponds to a region inside the unit circle [18]. This set of λ values was assigned by means of trial-and-error. The worst performance occurs for the highest setpoint values where the nonlinearities arise. The proposed approach provides a smoother and faster response than the typical PI with fixed parameters.To conclude, NG algorithm provides a robust procedure not only for clustering tasks, but also for feedforward nonlinear control using the gradient vectors obtained by the supervised version. These gradient vectors constitute the poles and zeros of the local transfer function of the plant. The computational complexity is linear regarding the number of samples, neurons, and variables because of the efficient implementation in batch procedure.



---

*Source: 3125073-2017-01-15.xml*