# MTN Optimal Tracking Control of SISO Nonlinear Time-Varying Discrete-Time Systems without Mechanism Models

Authors: Jiao-Jun Zhang; Hong-Sen Yan
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Nonlinear time-varying systems without mechanism models are common in application. They cannot be controlled directly by the traditional control methods based on precise mathematical models. Intelligent control is unsuitable for real-time control due to its computation complexity. For that sake, a multidimensional Taylor network (MTN) based output tracking control scheme, which consists of two MTNs, one as an identifier and the other as a controller, is proposed for SISO nonlinear time-varying discrete-time systems with no mechanism models. A MTN identifier is constructed to build the offline model of the system, and a set of initial parameters for online learning of the identifier is obtained. Then, an ideal output signal is selected relative to the given reference signal. Based on the system identification model, Pontryagin minimum principle is introduced to obtain the numerical solution of the optimal control law for the system relative to the given ideal output signal, with the corresponding optimal output taken as the desired output signal. A MTN controller is generated automatically to fit the numerical solution of the optimal control law using the conjugate gradient (CG) method, and a set of initial parameters for online learning of the controller is obtained. An adaptive back propagation (BP) algorithm is developed to adjust the parameters of the identifier and controller in real time, and the convergence for the proposed learning algorithm is verified. Simulation results show that the proposed scheme is valid.


---

## Body

## 1. Introduction

Nonlinear time-varying systems without mechanism models exist in practical engineering applications widely. However, it is difficult to obtain the precise mathematical model of a system due to the limitation of the modeling theory, the influences of its internal structure and parameter variations, and the external environment disturbances. In addition, the state variables are not easy to be determined, and it is inconvenient for state feedback control to be realized physically due to the practical and economic limitations of the measuring equipment in engineering practices. Output feedback control [1–6], which is of great theoretic and realistic significance, is to probe into the problem for nonlinear time-varying systems without mechanism models.Nonlinear autoregressive moving average with exogenous inputs (NARMAX) model describes an input-output relationship for a nonlinear dynamic system, by which the system output can be represented as a nonlinear functional expansion of its lagged inputs and outputs [7–9]. NARMAX has attracted considerable interest both in its theory and in applications [10–13], especially in the field of black-box nonlinear modeling. It is also referred to as time-varying NARMAX model in [14]. Neural network (NN) has the ability to approximate any continuous function with an arbitrary degree of accuracy over a compact set [15], and various kinds of NN have been used in system identification and control [7, 8, 16–24]. However, to satisfy the approximation requirement of a high-order uncertain system, both the numbers of the hidden layer neurons and the corresponding weight parameters needing online updating are large, which leads to the fact that the learning time tends to be unacceptably long and the real-time control is hardly realizable in practice. In addition, NN can only represent local dynamic characteristics when the state goes without too much change. However, the actual state change can be quite notable or even divergent. Thus NN cannot represent dynamic characteristics in the general sense. In fact, it cannot approximate to the unstable system as the artificial neurons suffer the limitations from the sign function, sigmoid function, and radial basis function, regardless of sample size or weight parameters. Furthermore, NN cannot approximate any continuous function to an arbitrarily degree of accuracy if the sign function or the sigmoid function is removed. From the point of view of the frequency characteristics of the input signals, a signal can be viewed as the superposition of value-high signal with the lowest frequency and value-low signals with high frequencies. Therefore, the value-low signals with high frequencies tend to be limited or restrained by the sign function, sigmoid function, or radial basis function of the artificial neuron, as a result of which the output of NN may fail to track the rapid change of its input. In terms of the neural network function approximation theorem, a three-layer neural network can approximate any nonlinear real-valued continuous function defined on a closed bounded subset. However, it cannot ensure that the nonlinear real-valued continuous function can be well approximated outside the closed bounded subset. Therefore, the general nonlinear dynamic system in the entire state space (i.e., not the bounded subset) is difficult to be approximated with an arbitrary degree of accuracy. As the polynomial function tends to be of infinity, the multidimensional Taylor network (MTN, whose idea was proposed by Hong-Sen Yan in 2010 and realization was done by Bo Zhou) is good at approximating or representing the general nonlinear dynamic system. It is suitable to be used as the identified dynamic model of the controlled plant, as it can represent polynomial dynamic system accurately, being simple with only states and inputs, and can be easily analyzed and solved for optimal control in terms of the minimum principle. However, NN only approximates the polynomial dynamic system instead of representing it. It is too complicated to be analyzed or solved for optimal control [25]. In addition, exponential function is contained in NN, which leads to the computational complexity and poor real-time control performance by a single chip microcomputer (SCM) and embedded system. That makes us resort to MTN, whereby only addition and multiplication are needed, its computation complexity being nearly that of the Taylor expansion of a single neuro in NN.MTN, first presented in [26], can reflect the dynamic characteristics of the system without knowing the order or other prior knowledge of the system. It approximates any nonlinear function with an arbitrary high accuracy, thus widely being applied in the study of time series prediction problems successfully [27–32]. The idea of MTN optimal control was proposed by Hong-Sen Yan in 2010 [33]. The optimal adjustment control of SISO nonlinear time-invariant systems has been achieved by introducing control input into MTN [34]. Asymptotic tracking and dynamic regulation of SISO nonlinear system based on discrete multidimensional Taylor network are considered in [35]. However, the system without mechanism model is not considered, nor are the time-varying characteristics of the nonlinear discrete-time system. MTN relies on the polynomial network for identification and control of the nonlinear system [10] in which the system considered is constant and its learning algorithm is based on the gradient descent algorithm with constant learning factor, which leads to its slow learning speed and convergence to local minima.Due to the uncertainties of the external environment and time-varying characteristics of a controlled plant [36, 37], the identifier and controller parameters need constant updating online in the process of its control, and the adjustments affect not only the control process but the robustness of the controller. Therefore, designing a desirable real-time self-tuning rule for the weight parameters of the identifier and controller is highly wanted. Backpropagation (BP) algorithm [38] is the most widely used learning algorithm in training multilayer neural network. However, it has such drawbacks as slow convergence speed and local optimal point. To raise the convergence speed, the improvement of BP algorithm [39–42] has been focused on, with certain desirable results achieved. However, in the improved learning algorithm, the learning rate and momentum factor are taken as constants in the interval (−1,1) randomly. There have emerged some other evolutionary algorithms developed to adjust the weight parameters for NN, such as genetic algorithm (GA) [43], particle swarm optimization (PSO) algorithm [44], the hybrid of them [45, 46], and the fuzzy logic approach [47], and so on. It has been shown that the coefficients should not remain fixed but should be changed adaptively throughout the entire training process so as to produce better training results, and it leads to the emergence of various schemes for adjusting the learning rate and momentum factor of BP algorithms adaptively [48, 49].In this paper, a MTN-based output tracking real-time control scheme, which consists of two MTNs, one as the identifier and the other as the controller, is proposed for SISO nonlinear time-varying discrete-time systems without mechanism models. A MTN identifier is developed for offline modeling of the system, and a set of initial parameters for online learning of the identifier is obtained. An ideal output signal is then selected for the given reference signal. Pontryagin minimum principle is employed to obtain the numerical solution of the optimal control law of the system relative to the ideal output signal, with the corresponding optimal output called the desired output signal. A MTN controller is generated automatically to fit the numerical solution of the optimal control law via the CG method, and a set of initial parameters for online learning of the controller is obtained. Based on the above, a novel adaptive BP algorithm for adjusting both the learning rate and the momentum factor in real time is designed to further enhance the learning speed of the identifier and controller. Finally, the convergence of the novel adaptive BP algorithm is analyzed. Simulation results show that the proposed scheme is valid.This paper is arranged as follows: in Section2, problem description; in Section 3, identifier design; in Section 4, automatic generation of the controller; in Section 5, selection of initial value of the online controller parameters; in Section 6, controller parameters self-tuning; The algorithm steps for MTN optimal control scheme are summarized in Section 7; simulation study is mentioned in Section 8; conclusion is in Section 9.


## 2. Problem Statement

Consider the following unknown SISO nonlinear time-varying discrete-time system described by the input-output difference equation:(1)yk+1=fyk,yk-1,…,yk-ny,uk,uk-1,…,uk-nu,k,where f(·) is an unknown nonlinear scale function, y(k) and u(k) are the output and input of the system, and ny and nu are the corresponding maximum delays.The goals of the present study are as follows: (a) to design an offline identifier to build the system model based on the input-output data pairs{u(k),y(k+1)}; (b) to design such real-time controller that allows the output y(k) of the system (1) to track the given reference signal r(k) as closely as possible.The block diagram of control system (1) is shown in Figure 1. For clarity, the multidimensional Taylor network identifier and controller are abbreviated as MTNI and MTNC, as shown in Figure 1.Figure 1
Block diagram of control system (1).


## 3. System Identification

It is known from [26] that MTN provides a good nonlinear function approximation approach, and f(·) in the system (1) can be approximated with arbitrary precision by MTN, using an appropriate learning algorithm. Let f^(·) be the mapping relationship, and we obtain the MTN model (MTNI) of the system (1) as follows:(2)y^k+1=f^k^1yk,k^2yk-1,…,k^ny+1yk-ny,l^1uk,l^2uk-1,…,l^nu+1uk-nu,w^kwhere y^(k+1) is the output of MTNI, w^(k) is the weight coefficient vector of MTNI, and l^i and k^j are positive constants, i=1,2,…,nu+1, j=1,2,…,ny+1.For convenience and without loss of generality, sett^=ny+nu+2 and we get(3)z^k=z^1k,…,z^ny+1k,z^ny+2k,…,z^t^-1k,z^t^kT=k^1yk,…,k^ny+1yk-ny,l^2uk-1,…,l^nu+1uk-nu,l^1ukT.Setting the weight coefficient vector of MTNI asw^(k)=[w^1(k),…,w^N^(t^,m^)(k)]T allows us to rewrite the identification model (2) as(4)y^k+1=∑p^=1N^t^,m^w^p^k∏q^=1t^z^q^λ^p^,q^k,where N^(t^,m^) represents the total number of product items of the t^-ary function f^(·) expanded into the approximate polynomial with m^ powers, w^p^(k) is the weight coefficient of the p^-th product item in the formula, λ^(p^,q^) denotes the power of the variable z^q^(k) in the p^-th product item, and ∑q^=1t^λ^(p^,q^)≤m^, where p^=1,2,…,N^(t^,m^).The diagram of MTNI is shown in Figure2.Figure 2
The diagram of MTNI.To calculateN^(t^,m^) and λ^(p^,q^), the product items in (4) are rearranged as illustrated in Figure 3, i.e., storing the product items of the expansion according to their powers, respectively. We use the symbol (i^,j^) to denote the i^-th rectangle in which the product items with j^-th power are stored and store the product items with j^-th power which are got by adding 1 on the power of the i^-th element z^i^(k) from the i^-th rectangle to the t^-th rectangle with j^-1-th power into (i^,j^), and so on, until storing the product items with m^-th power which are obtained by adding 1 on the power of the t^-th element z^t^(k) in (t^,m^-1) into (t^,m^), where i^=1,2,…,t^ and j^=2,3,…,m^.Figure 3
The new arrangement of product items for MTNI.The calculation ofN^(t^,m^) and λ^(p^,q^) goes as follows.LetP^(i^,j^) represent the number of product items in (i^,j^), and we get(5)N^t^,m^=∑j^=1m^∑i^=1t^P^i^,j^+1where(6)P^i^,j^=∑s^=i^t^P^s^,j^-1,j^=2,3,…,m^,1,j^=1.Suppose that, in (4), from the 2-th item, the p^-th product item corresponds to the r^-th product item in (i^,j^) of Figure 3. For clarity, the power of the element z^q^(k) is termed as λ^i^,j^(r^,q^), and Q^j^(a^,b^) represent the number of product items with the j^-th power from the a^-th to the b^-th rectangle. From Figure 3, it is known that(7)λ^i^,j^r^,q^=λ^i^+d^,j^-1r^-Q^j^-1i^,i^+d^-1,q^+1,q^=i^,λ^i^+d^,j^-1r^-Q^j^-1i^,i^+d^-1,q^,q^≠i^,j^=2,3,…,m^Q^j^a^,b^=∑s^=a^b^P^s^,j^,a^≤b^,0,a^>b^,j^=2,3,…,m^.The initial values are set as follows:(8)λ^i^,j^r^,q^=1,q^=i^,0,q^≠i^,j^=1,Q^j^a^,b^=b^-a^+1,a^≤b^,0,a^>b^,j^=1,where Q^j^-1(i^,i^+d^-1)<r^≤Q^j^-1(i^,i^+d^),d^=0,1,…,t^-i^ and y=y(1),y(2),…,y(N)T.Based on the above, Figure4 gives the diagram of system identification.Figure 4
The diagram of system identification.


### 3.1. Offline System Identification

The identification error can be defined as(9)e⌢k+1=y^k+1-yk+1,where k=0,1,…,N-1.The corresponding mean square error is(10)E⌢=1N∑k=0N-1e⌢2k+1.Substituting (4) and (9) into (10) yields(11)E⌢=1N∑k=0N-1∑p^=1N^t^,m^w^p^∏q^=1t^z^q^λ^p^,q^k-yk+12.Assume(12)α⌢k=∏q^=1t^z^q^λ^1,q^k,∏q^=1t^z^q^λ^2,q^k,…,∏q^=1t^z^q^λ^N^t^,m^,q^kTA⌢=α⌢0,α⌢1,…,α⌢N-1Equation (11) can be rewritten into(13)E⌢=1N∑k=0N-1w^Tα⌢k-yk+12=1Nw^T∑k=0N-1α⌢kα⌢Tkw^-2N∑k=0N-1yk+1α⌢Tkw^+1N∑k=0N-1y2k+1.Calculate the partial derivative ofE⌢ with respect to the weight coefficient vector w^, i.e.,(14)∂E⌢∂w^=2N∑k=0N-1α⌢kα⌢Tkw^-2N∑k=0N-1yk+1α⌢k.Settingg⌢=∂E⌢/∂w^, Q⌢=2/N∑k=0N-1α⌢kα⌢Tk and b⌢=-2/N∑k=0N-1y(k+1)α⌢(k) gives(15)g⌢=Q⌢w^+b⌢.Lettingy=y(1),y(2),…,y(N)T, the vector form of Q⌢, b⌢ and g⌢ can be rewritten as(16)Q⌢=2NA⌢A⌢T,b⌢=-2NA⌢y,g⌢=2NA⌢A⌢Tw^-y.To obtain a precise model of system (1), the weight coefficient vector w^ should be updated over and over again by observation of the input-output data pairs {u(k),y(k+1)}. A number of classical weight update laws have been proposed in the literatures, such as least squares algorithms, various gradient-type algorithms [50, 51], least-mean-square (LMS) algorithm [52], etc. The gradient method is commonly adopted for parameter adjustment; that is, w^ can be updated once in the negative gradient direction after each offline learning. Let w^τ⌢ represent the value of w^ after the τ⌢ th training, and we obtain w^τ⌢+1=w^τ⌢-μ⌢τ⌢g⌢τ⌢, where μ⌢τ⌢=g⌢τ⌢Tg⌢τ⌢/g⌢τ⌢TQ⌢g⌢τ⌢. However, as the gradient path to the minimum point is zigzag, the search direction remains vertical to the last. Fortunately, the problem can be solved effectively by employing the CG method, whereby the weight can be updated as follows:(17)w^τ⌢+1=w^τ⌢+μ⌢τ⌢p⌢τ⌢.where p⌢τ⌢=-g⌢τ⌢+a⌢τ⌢-1p⌢τ⌢-1, μ⌢τ⌢=g⌢τ⌢Tg⌢τ⌢/p⌢τ⌢TQ⌢p⌢τ⌢, a⌢τ⌢-1=g⌢τ⌢TQ⌢p⌢τ⌢-1/p⌢τ⌢-1TQ⌢p⌢τ⌢-1, and the initial value is p⌢1=-g⌢1.The block diagram of offline learning of MTNI (4) is shown in Figure 5.Figure 5
Block diagram of offline learning of MTNI.


### 3.2. Real-Time Learning for the Weights of MTNI

For an unknown system, system identification is the mathematical modeling process by observation of the input-output data pairs. Nonlinear time-varying system identification based on MTNI is to take the connection weight coefficients of MTNI as time-varying parameters to be estimated and trained online by suitable learning algorithm, with the same outputs of the plant and the model for the same set of inputs. The weight coefficients need to be adjusted online for desirable real-time identification effect.Set the performance function for MTNI as(18)E^w^k=12e^2w^k,where e^(w^(k))=e^(k+1),  e^(k+1)=y(k+1)-y^(k+1) represents the identification error at time k+1, and k=0,1,….To obtain a better identification effect for the unknown nonlinear time-varying system, the weight coefficients of MTNI should be adjusted adaptively throughout the entire training process. A novel adaptive BP algorithm for adjusting both the learning rate and momentum factor adaptively [53–55] is proposed, and the weight parameters are updated, i.e.,(19)w^k+1=w^k+Δw^k,where(20)Δw^k=-η^k∇E^w^k+α^kΔw^k-1,(21)η^k=min⁡η^k-11+η^0cos⁡θ^k,η^max,(22)α^k=β^kη^2k,(23)β^k=β^0-∇E^w^k·Δw^k-1Δw^k-12,if∇E^w^k·Δw^k-1<0,0,else,where ∇E^(w^(k)) is the partial derivative of E^(w^(k)) with respect to w^(k); η^0, η^max, β^0 are constants, and 0≤η^0≤1, 0<η^max<1, 0≤β^0<1; θ^(k) is the angle between the current gradient -∇E^(w^(k)) and previous update Δw^(k-1), given by cos⁡θ^(k)=-∇E^(w^(k))·Δw^(k-1)/∇E^(w^(k))Δw^(k-1).Theorem 1.
For any given set of weight coefficient vectorw^(0), if w^(k) is generated by the learning rules from (19) to (23), there exists ∇E^(w^(k))·Δw^(k)≤0.Proof.
For the first case, i.e.,∇E^(w^(k))·Δw^(k-1)<0, we have(24)
Δ
w
^
k
=
-
η
^
k
∇
E
^
w
^
k
+
α
^
k
Δ
w
^
k
-
1
=
η
^
k
-
∇
E
^
w
^
k
+
β
^
0
η
^
k
-
∇
E
^
w
^
k
·
Δ
w
^
k
-
1
Δ
w
^
k
-
1
2
Δ
w
^
k
-
1
=
η
^
k
-
∇
E
^
w
^
k
+
β
^
0
η
^
k
∇
E
^
w
^
k
Δ
w
^
k
-
1
cos
⁡
θ
^
k
Δ
w
^
k
-
1
2
Δ
w
^
k
-
1
=
η
^
k
-
∇
E
^
w
^
k
+
β
^
0
η
^
k
cos
⁡
θ
^
k
∇
E
^
w
^
k
Δ
w
^
k
-
1
Δ
w
^
k
-
1
,
(25)
∇
E
^
w
^
k
·
Δ
w
^
k
=
∇
E
^
T
w
^
k
Δ
w
^
k
=
η
^
k
-
∇
E
^
w
^
k
2
+
β
^
0
η
^
k
cos
⁡
θ
^
k
∇
E
^
w
^
k
Δ
w
^
k
-
1
∇
E
^
T
w
^
k
Δ
w
^
k
-
1
=
-
η
^
k
∇
E
^
w
^
k
2
1
+
β
^
0
η
^
k
cos
2
⁡
θ
^
k
As a matter of fact, there existη^(k)≥0 and 1+β^0η^kcos2⁡θ^k>0 when 0≤β^0<1 and 0≤η^0≤1, as a result of which ∇E^(w^(k))·Δw^(k)≤0 holds.
For the second case, i.e.,∇E^(w^(k))·Δw^(k-1)≥0, we get β^(k)=0, α^(k)=0 and Δw^(k)=-η^(k)∇E^(w^(k)), thus,(26)∇E^w^k·Δw^k=∇E^Tw^kΔw^k=-η^k∇E^w^k2≤0.
As revealed by the above two cases,∇E^(w^(k))·Δw^(k)≤0 holds.
That completes the proof of Theorem1.


## 3.1. Offline System Identification

The identification error can be defined as(9)e⌢k+1=y^k+1-yk+1,where k=0,1,…,N-1.The corresponding mean square error is(10)E⌢=1N∑k=0N-1e⌢2k+1.Substituting (4) and (9) into (10) yields(11)E⌢=1N∑k=0N-1∑p^=1N^t^,m^w^p^∏q^=1t^z^q^λ^p^,q^k-yk+12.Assume(12)α⌢k=∏q^=1t^z^q^λ^1,q^k,∏q^=1t^z^q^λ^2,q^k,…,∏q^=1t^z^q^λ^N^t^,m^,q^kTA⌢=α⌢0,α⌢1,…,α⌢N-1Equation (11) can be rewritten into(13)E⌢=1N∑k=0N-1w^Tα⌢k-yk+12=1Nw^T∑k=0N-1α⌢kα⌢Tkw^-2N∑k=0N-1yk+1α⌢Tkw^+1N∑k=0N-1y2k+1.Calculate the partial derivative ofE⌢ with respect to the weight coefficient vector w^, i.e.,(14)∂E⌢∂w^=2N∑k=0N-1α⌢kα⌢Tkw^-2N∑k=0N-1yk+1α⌢k.Settingg⌢=∂E⌢/∂w^, Q⌢=2/N∑k=0N-1α⌢kα⌢Tk and b⌢=-2/N∑k=0N-1y(k+1)α⌢(k) gives(15)g⌢=Q⌢w^+b⌢.Lettingy=y(1),y(2),…,y(N)T, the vector form of Q⌢, b⌢ and g⌢ can be rewritten as(16)Q⌢=2NA⌢A⌢T,b⌢=-2NA⌢y,g⌢=2NA⌢A⌢Tw^-y.To obtain a precise model of system (1), the weight coefficient vector w^ should be updated over and over again by observation of the input-output data pairs {u(k),y(k+1)}. A number of classical weight update laws have been proposed in the literatures, such as least squares algorithms, various gradient-type algorithms [50, 51], least-mean-square (LMS) algorithm [52], etc. The gradient method is commonly adopted for parameter adjustment; that is, w^ can be updated once in the negative gradient direction after each offline learning. Let w^τ⌢ represent the value of w^ after the τ⌢ th training, and we obtain w^τ⌢+1=w^τ⌢-μ⌢τ⌢g⌢τ⌢, where μ⌢τ⌢=g⌢τ⌢Tg⌢τ⌢/g⌢τ⌢TQ⌢g⌢τ⌢. However, as the gradient path to the minimum point is zigzag, the search direction remains vertical to the last. Fortunately, the problem can be solved effectively by employing the CG method, whereby the weight can be updated as follows:(17)w^τ⌢+1=w^τ⌢+μ⌢τ⌢p⌢τ⌢.where p⌢τ⌢=-g⌢τ⌢+a⌢τ⌢-1p⌢τ⌢-1, μ⌢τ⌢=g⌢τ⌢Tg⌢τ⌢/p⌢τ⌢TQ⌢p⌢τ⌢, a⌢τ⌢-1=g⌢τ⌢TQ⌢p⌢τ⌢-1/p⌢τ⌢-1TQ⌢p⌢τ⌢-1, and the initial value is p⌢1=-g⌢1.The block diagram of offline learning of MTNI (4) is shown in Figure 5.Figure 5
Block diagram of offline learning of MTNI.


## 3.2. Real-Time Learning for the Weights of MTNI

For an unknown system, system identification is the mathematical modeling process by observation of the input-output data pairs. Nonlinear time-varying system identification based on MTNI is to take the connection weight coefficients of MTNI as time-varying parameters to be estimated and trained online by suitable learning algorithm, with the same outputs of the plant and the model for the same set of inputs. The weight coefficients need to be adjusted online for desirable real-time identification effect.Set the performance function for MTNI as(18)E^w^k=12e^2w^k,where e^(w^(k))=e^(k+1),  e^(k+1)=y(k+1)-y^(k+1) represents the identification error at time k+1, and k=0,1,….To obtain a better identification effect for the unknown nonlinear time-varying system, the weight coefficients of MTNI should be adjusted adaptively throughout the entire training process. A novel adaptive BP algorithm for adjusting both the learning rate and momentum factor adaptively [53–55] is proposed, and the weight parameters are updated, i.e.,(19)w^k+1=w^k+Δw^k,where(20)Δw^k=-η^k∇E^w^k+α^kΔw^k-1,(21)η^k=min⁡η^k-11+η^0cos⁡θ^k,η^max,(22)α^k=β^kη^2k,(23)β^k=β^0-∇E^w^k·Δw^k-1Δw^k-12,if∇E^w^k·Δw^k-1<0,0,else,where ∇E^(w^(k)) is the partial derivative of E^(w^(k)) with respect to w^(k); η^0, η^max, β^0 are constants, and 0≤η^0≤1, 0<η^max<1, 0≤β^0<1; θ^(k) is the angle between the current gradient -∇E^(w^(k)) and previous update Δw^(k-1), given by cos⁡θ^(k)=-∇E^(w^(k))·Δw^(k-1)/∇E^(w^(k))Δw^(k-1).Theorem 1.
For any given set of weight coefficient vectorw^(0), if w^(k) is generated by the learning rules from (19) to (23), there exists ∇E^(w^(k))·Δw^(k)≤0.Proof.
For the first case, i.e.,∇E^(w^(k))·Δw^(k-1)<0, we have(24)
Δ
w
^
k
=
-
η
^
k
∇
E
^
w
^
k
+
α
^
k
Δ
w
^
k
-
1
=
η
^
k
-
∇
E
^
w
^
k
+
β
^
0
η
^
k
-
∇
E
^
w
^
k
·
Δ
w
^
k
-
1
Δ
w
^
k
-
1
2
Δ
w
^
k
-
1
=
η
^
k
-
∇
E
^
w
^
k
+
β
^
0
η
^
k
∇
E
^
w
^
k
Δ
w
^
k
-
1
cos
⁡
θ
^
k
Δ
w
^
k
-
1
2
Δ
w
^
k
-
1
=
η
^
k
-
∇
E
^
w
^
k
+
β
^
0
η
^
k
cos
⁡
θ
^
k
∇
E
^
w
^
k
Δ
w
^
k
-
1
Δ
w
^
k
-
1
,
(25)
∇
E
^
w
^
k
·
Δ
w
^
k
=
∇
E
^
T
w
^
k
Δ
w
^
k
=
η
^
k
-
∇
E
^
w
^
k
2
+
β
^
0
η
^
k
cos
⁡
θ
^
k
∇
E
^
w
^
k
Δ
w
^
k
-
1
∇
E
^
T
w
^
k
Δ
w
^
k
-
1
=
-
η
^
k
∇
E
^
w
^
k
2
1
+
β
^
0
η
^
k
cos
2
⁡
θ
^
k
As a matter of fact, there existη^(k)≥0 and 1+β^0η^kcos2⁡θ^k>0 when 0≤β^0<1 and 0≤η^0≤1, as a result of which ∇E^(w^(k))·Δw^(k)≤0 holds.
For the second case, i.e.,∇E^(w^(k))·Δw^(k-1)≥0, we get β^(k)=0, α^(k)=0 and Δw^(k)=-η^(k)∇E^(w^(k)), thus,(26)∇E^w^k·Δw^k=∇E^Tw^kΔw^k=-η^k∇E^w^k2≤0.
As revealed by the above two cases,∇E^(w^(k))·Δw^(k)≤0 holds.
That completes the proof of Theorem1.


## 4. Controller

The control objective is to find a control inputu(k) that enables the system output to track in real time the given reference signal r(k) as closely as possible in real time. In this section, we consider the controller MTNC generated automatically as follows:(27)uk=hk1uk-1,…,kn-uuk-n-u,l1ek-1,…,lneek-ne,k,where u(k) is the output of MTNC, i.e., the input of the system (1), e(k)=r(k)-y(k) is the tracking error at time k, n-u and ne are the maximum delays of the output and input of MTNC, and li and kj are positive constants, i=1,2,…,ne, j=1,2,…,n-u.For convenience, without loss of generality, denotet=n-u+ne, and we have(28)zk=z1k,…,zn-uk,zn-u+1k,…,ztkT=k1uk-1,…,kn-uuk-n-u,l1ek-1,…,lneek-neT.Known from [26], there exists a group of weight coefficient vectors w(k)=[w1(k),…,wN(t,m)(k)]T, and thus, the input u(k) can be rewritten as(29)uk=∑p=1Nt,mwpk∏q=1tzqλp,qk,where N(t,m) represents the total number of product items for the t-ary function h(·) expanded into the approximate polynomial with m powers, wp(k) denotes the weight coefficient of the p-th product item, λ(p,q) is the power of the variable zq(k) in the p-th product item, and ∑q=1tλ(p,q)≤m, where p=1,2,…,N(t,m).The diagram of MTNC is shown in Figure6.Figure 6
The diagram of MTNC.To calculateN(t,m) and λ(p,q), the product items in (29) are rearranged as illustrated in Figure 7, i.e., storing the product items of the expansion according to their powers, respectively. We use the symbol (i,j) to denote the i-th rectangle in which the product items with j-th power are stored and store the product items with j-th power which are got by adding 1 on the power of the i-th element zi(k) from the i-th rectangle to the t-th rectangle with j-1-th power into (i,j), and so on, until storing the product items with m-th power which are obtained by adding 1 on the power of the t-th element zt(k) in (t,m-1) into (t,m), where i=1,2,…,t and j=2,3,…,m.Figure 7
The new arrangement of the product items for MTNC.The calculation ofN(t,m) and λ(p,q) goes as follows.LetP(i,j) represent the number of product items in (i,j), and we get(30)Nt,m=∑j=1m∑i=1tPi,j+1,where(31)Pi,j=∑s=itPs,j-1,j=2,3,…,m,1,j=1.Suppose that, in (29), from the 2-th item, the p-th product item is according to the r-th product item in (i,j) of Figure 7. For clarity, the power of the element zq(k) is termed as λi,j(r,q), and Qj(a,b) represents the number of product items with the j-th power from the a-th to the b-th rectangle. From Figure 7, it is known that(32)λi,jr,q=λi+d,j-1r-Qj-1i,i+d-1,q+1,q=i,λi+d,j-1r-Qj-1i,i+d-1,q,q≠i,Qja,b=∑s=abPs,j,a≤b,0,a>b.The initial values are set as follows:(33)λi,jr,q=1,q=i,0,q≠i,j=1,Qja,b=b-a+1,a≤b,0,a>b,j=1,where Qj-1(i,i+d-1)<r≤Qj-1(i,i+d), d=0,1,…,t-i and i=1,2,…,t.


## 5. Initial Weight Values of MTNC

The convergence speed is influenced by the selection of the initial weight values of the controller. However, random choice of network parameters is the most common practice in network training. To enhance the convergence speed and avoid falling into local minimum, two steps are introduced here for selection of the initial weight values. The first is to transform the offline model (4) of the system into an extended state space description form through variable substitution, select a group of ideal output signal yid(k) relative to the given reference signal r(k), and employ Pontryagin’s minimum principle to obtain the numerical solution of the optimal control law u∗(k) of the system relative to the ideal output signal yid(k), with the corresponding optimal output called desired output signal yop(k). In the second step, a set of parameter values is given randomly in the interval (-1,1), and the CG method is applied for MTNC offline training to approximate the optimal control law u∗(k). A set of weight values are then obtained as the initial values for online training MTNC, where k=0,1,…,N-1. The specific steps go as follows.


### 5.1. Optimal Control Law

Based on the identification model (4) obtained offline for the system (1), y(k) can be substituted by y^(k). For convenience, set(34)xk=x1k,x2k,…,xny+1k,xny+2k,…,xt^-1kT=k^1y^k,k^2y^k-1,…,k^ny+1y^k-ny,l^2uk-1,…,l^nu+1uk-nuT,then, we obtain the extended state space description form with the following variable substitution:(35)x1k+1=k^1∑p^=1N^t^,m^w^p^∏q^=1t^-1xq^λ^p^,q^kl^1λ^p^,t^uλ^p^,t^kxjk+1=k^jk^j-1xj-1k,2≤j≤ny+1xny+2k+1=l^2ukxjk+1=l^j-nyl^j-1-nyxj-1k,ny+3≤j≤t^-1(36)y^k=1k^1x1k.Consider the following optimal control problem [56]:(37)minukE~=12∑k=0N-1yidk-y^k2,where y^(k) satisfies such the constraint conditions as (35) and (36).Introduce the Hamiltonian equation:(38)Hk=Hxk,uk,λk+1,k=12yidk-1k^1x1k2+λ1k+1k^1∑p^=1N^t^,m^w^p^∏q^=1t^-1xq^λ^p^,q^kl^1λ^p^,t^uλ^p^,t^k+∑j=2ny+1λjk+1k^jk^j-1xj-1k+λny+2k+1l^2uk+∑j=ny+3t^-1λjk+1l^j-nyl^j-1-nyxj-1kwhere x(k) and λ(k) satisfy the following conditions:(39)xk+1=∂Hk∂λk+1,(40)λk=∂Hk∂xk.If the control vector is constrained, Hamiltonian function takes the extreme value on the optimal control sequence by the minimum principle; i.e., take extreme value on the extreme valuesx∗(k) of the optimal trajectory and the optimal control law u∗(k), that is,(41)Hx∗k,u∗k,λk+1,k=minuk∈ΩHx∗k,uk,λk+1,k,where Ω is a bounded closed set.If the control vector is not constrained, Hamiltonian function takes extreme value from the whole control spaceR, the extreme condition being(42)∂Hk∂uk=0,and u(k)∈R.The given series of control sequenceu(k)(k=0,1,…,N-1) can be improved by repeated iteration in the direction that makes the gradient of Hamiltonian function H(k) decrease, until the necessary condition (42) is satisfied. Then we obtain the numerical solution of the optimal control law u∗(k), where k=0,1,…,N-1. For convenience, let U=(u(0),u(1),…,u(N-1))T, and the calculation steps are as follows.Algorithm 2.
⁡
Step  1. Set any given series of control sequence uM(k)=u0(k), where M is the number of iterations, and the initial value is set as M=0, and k=0,1,…,N-1.
Step  2. Solve the state variable xM(k) sequentially by formula (35) based on UM and the initial condition x(0), where k=1,2,…,N.
Step  3. Calculate gM(k)=∂Hk/∂uku(k)=uM(k), which is the gradient of H(k) with respect to u(k) in the control sequence uM(k), and set GM=(gM(0),gM(1),…,gM(N-1))T, where k=0,1,…,N-1.
Step  4. Calculate GM. If GM<ε, stop, or else, revise the control vector: UM+1=UM-σGM, i.e., uM+1(k)=uM(k)-σgM(k), where ε is a given value, σ is a fixed step size, and k=0,1,…,N-1.
Step  5. Let M=M+1, and return to Step 2.


### 5.2. Initial Weight Values of MTNC

MTNC is generated automatically to approximate the numerical solutionu∗(k) of the optimal control law for the offline model (4) of the system (1) relative to the ideal output signal yid(k). The weight coefficients of MTNC are obtainable by offline learning, and the block diagram of offline learning for MTNC (29) is shown in Figure 8.Figure 8
Block diagram of offline learning of MTNC.The initial weight values of MTNC are secured offline by the CG method as follows.Define the appropriate error as(43)e-k=u∗k-uk,where u(k) is the output of the controller MTNC at time k.The corresponding mean square error is(44)E-=12∑k=0N-1e-2k.Substituting (43) and (29) into (44) yields(45)E-=12∑k=0N-1u∗k-∑p=1Nt,mwp∏q=1tzqλp,qk2.Set(46)αk=∏q=1tzqλ1,qk,∏q=1tzqλ2,qk,…,∏q=1tzqλNt,m,qkT,A=α0,α1,…,αN-1,and formula (45) can be rewritten as(47)E-=12∑k=0N-1u∗k-wTαk2=12wT∑k=0N-1αkαTkw-∑k=0N-1u∗kαTkw+12∑k=0N-1u∗2k.Calculate the partial derivative ofE- with respect to the weight coefficient vector w:(48)∂E-∂w=∑k=0N-1αkαTkw-∑k=0N-1u∗kαk.Letg=∂E-/∂w, Q=∑k=0N-1(α(k)αT(k)), and b=-∑k=0N-1u∗(k)α(k), and we have(49)g=Qw+b.Settingu∗=(u∗(0),u∗(1),…,u∗(N-1))T enables us to get the vector form of Q, b and g as follows:(50)Q=AAT,b=-Au∗,g=AATw-u∗.For the given numerical solutionu∗(k) of the optimal control, the weight coefficient vector w can be updated once in the negative gradient direction after each learning. Let wτ represent the value of w after the τth iterative training; then we have wτ+1=wτ-μτgτ, where μτ=gτTgτ/gτTQgτ. However, as the gradient path to the minimum point is zigzag, the search direction remains vertical to the last. Fortunately, the problem can be solved effectively by employing the CG method, whereby the weight can be updated as follows:(51)wτ+1=wτ+μτpτ,where pτ=-gτ+aτ-1pτ-1, μτ=gτTgτ/pτTQpτ, aτ-1=gτTQpτ-1/pτ-1TQpτ-1, and the initial value is defined as p1=-g1.


## 5.1. Optimal Control Law

Based on the identification model (4) obtained offline for the system (1), y(k) can be substituted by y^(k). For convenience, set(34)xk=x1k,x2k,…,xny+1k,xny+2k,…,xt^-1kT=k^1y^k,k^2y^k-1,…,k^ny+1y^k-ny,l^2uk-1,…,l^nu+1uk-nuT,then, we obtain the extended state space description form with the following variable substitution:(35)x1k+1=k^1∑p^=1N^t^,m^w^p^∏q^=1t^-1xq^λ^p^,q^kl^1λ^p^,t^uλ^p^,t^kxjk+1=k^jk^j-1xj-1k,2≤j≤ny+1xny+2k+1=l^2ukxjk+1=l^j-nyl^j-1-nyxj-1k,ny+3≤j≤t^-1(36)y^k=1k^1x1k.Consider the following optimal control problem [56]:(37)minukE~=12∑k=0N-1yidk-y^k2,where y^(k) satisfies such the constraint conditions as (35) and (36).Introduce the Hamiltonian equation:(38)Hk=Hxk,uk,λk+1,k=12yidk-1k^1x1k2+λ1k+1k^1∑p^=1N^t^,m^w^p^∏q^=1t^-1xq^λ^p^,q^kl^1λ^p^,t^uλ^p^,t^k+∑j=2ny+1λjk+1k^jk^j-1xj-1k+λny+2k+1l^2uk+∑j=ny+3t^-1λjk+1l^j-nyl^j-1-nyxj-1kwhere x(k) and λ(k) satisfy the following conditions:(39)xk+1=∂Hk∂λk+1,(40)λk=∂Hk∂xk.If the control vector is constrained, Hamiltonian function takes the extreme value on the optimal control sequence by the minimum principle; i.e., take extreme value on the extreme valuesx∗(k) of the optimal trajectory and the optimal control law u∗(k), that is,(41)Hx∗k,u∗k,λk+1,k=minuk∈ΩHx∗k,uk,λk+1,k,where Ω is a bounded closed set.If the control vector is not constrained, Hamiltonian function takes extreme value from the whole control spaceR, the extreme condition being(42)∂Hk∂uk=0,and u(k)∈R.The given series of control sequenceu(k)(k=0,1,…,N-1) can be improved by repeated iteration in the direction that makes the gradient of Hamiltonian function H(k) decrease, until the necessary condition (42) is satisfied. Then we obtain the numerical solution of the optimal control law u∗(k), where k=0,1,…,N-1. For convenience, let U=(u(0),u(1),…,u(N-1))T, and the calculation steps are as follows.Algorithm 2.
⁡
Step  1. Set any given series of control sequence uM(k)=u0(k), where M is the number of iterations, and the initial value is set as M=0, and k=0,1,…,N-1.
Step  2. Solve the state variable xM(k) sequentially by formula (35) based on UM and the initial condition x(0), where k=1,2,…,N.
Step  3. Calculate gM(k)=∂Hk/∂uku(k)=uM(k), which is the gradient of H(k) with respect to u(k) in the control sequence uM(k), and set GM=(gM(0),gM(1),…,gM(N-1))T, where k=0,1,…,N-1.
Step  4. Calculate GM. If GM<ε, stop, or else, revise the control vector: UM+1=UM-σGM, i.e., uM+1(k)=uM(k)-σgM(k), where ε is a given value, σ is a fixed step size, and k=0,1,…,N-1.
Step  5. Let M=M+1, and return to Step 2.


## 5.2. Initial Weight Values of MTNC

MTNC is generated automatically to approximate the numerical solutionu∗(k) of the optimal control law for the offline model (4) of the system (1) relative to the ideal output signal yid(k). The weight coefficients of MTNC are obtainable by offline learning, and the block diagram of offline learning for MTNC (29) is shown in Figure 8.Figure 8
Block diagram of offline learning of MTNC.The initial weight values of MTNC are secured offline by the CG method as follows.Define the appropriate error as(43)e-k=u∗k-uk,where u(k) is the output of the controller MTNC at time k.The corresponding mean square error is(44)E-=12∑k=0N-1e-2k.Substituting (43) and (29) into (44) yields(45)E-=12∑k=0N-1u∗k-∑p=1Nt,mwp∏q=1tzqλp,qk2.Set(46)αk=∏q=1tzqλ1,qk,∏q=1tzqλ2,qk,…,∏q=1tzqλNt,m,qkT,A=α0,α1,…,αN-1,and formula (45) can be rewritten as(47)E-=12∑k=0N-1u∗k-wTαk2=12wT∑k=0N-1αkαTkw-∑k=0N-1u∗kαTkw+12∑k=0N-1u∗2k.Calculate the partial derivative ofE- with respect to the weight coefficient vector w:(48)∂E-∂w=∑k=0N-1αkαTkw-∑k=0N-1u∗kαk.Letg=∂E-/∂w, Q=∑k=0N-1(α(k)αT(k)), and b=-∑k=0N-1u∗(k)α(k), and we have(49)g=Qw+b.Settingu∗=(u∗(0),u∗(1),…,u∗(N-1))T enables us to get the vector form of Q, b and g as follows:(50)Q=AAT,b=-Au∗,g=AATw-u∗.For the given numerical solutionu∗(k) of the optimal control, the weight coefficient vector w can be updated once in the negative gradient direction after each learning. Let wτ represent the value of w after the τth iterative training; then we have wτ+1=wτ-μτgτ, where μτ=gτTgτ/gτTQgτ. However, as the gradient path to the minimum point is zigzag, the search direction remains vertical to the last. Fortunately, the problem can be solved effectively by employing the CG method, whereby the weight can be updated as follows:(51)wτ+1=wτ+μτpτ,where pτ=-gτ+aτ-1pτ-1, μτ=gτTgτ/pτTQpτ, aτ-1=gτTQpτ-1/pτ-1TQpτ-1, and the initial value is defined as p1=-g1.


## 6. Real-Time Learning for Weights of MTNC

Due to the real-time modeling error for the unknown nonlinear time-varying system and the uncertainties existing in practical applications, the controller with fixed weight coefficients cannot ensure the lasting robust performance of the system. Therefore, it is required that the controller be capable of adjusting automatically for real-time control. Similar to the real-time learning algorithm for MTNI, a novel adaptive BP algorithm is proposed here for MTNC real-time training [53–55].The performance function is defined as(52)Ewk=12e2wk,where e(w(k))=e(k+1), e(k+1)=r(k+1)-y(k+1) represents the practical tracking error at time k+1.The weight coefficients of MTNC are updated according to (53)-(57):(53)wk+1=wk+Δwk,where(54)Δwk=-ηk∇Ewk+αkΔwk-1,(55)ηk=min⁡ηk-11+η0cos⁡θk,ηmax,(56)αk=βkη2k,(57)βk=β0-∇Ewk·Δwk-1Δwk-12,if∇Ewk·Δwk-1<0,0,else,∇E(w(k)) is the partial derivative of E(w(k)) with respect to w(k); η0, ηmax, β0 are constants, and 0≤η0≤1, 0<ηmax<1, 0≤β0<1; θ(k) is the angle between the current gradient -∇E(w(k)) and the previous update Δw(k-1), given by cos⁡θ(k)=-∇E(w(k))·Δw(k-1)/∇E(w(k))Δw(k-1).From (54) and (57), we get(58)∇Ewk=ewkdewkdwk=-ewkdyk+1dwk.As the system considered is unknown, the outputy(k) of the actual system can be replaced with the output y^(k) of the identification model (4), that is,(59)yk+1≈y^k+1=∑p^=1N^t^,m^w^p^k∏q^=1t^z^q^λ^p^,q^k,where t^, m^, N^(t^,m^), p^, q^, λ^(p^,q^), w^p^(k) and z^q^(k) have the same meanings as in (4), and k=0,1,2,….Thus,(60)dyk+1dwk≈dy^k+1dwk,(61)dy^k+1dwk=l^1∂y^k+1∂z^t^kdukdwk.In (61), the first part on the right side of the equation can be calculated using the real-time identification model (4), and the second term can be done as follows:(62)dukdwk=αkwhere α(k) is mentioned before.Theorem 3.
For any given set of weight vectorw0 of MTNC used to approximate the optimal control law u∗(k) offline, by taking the weight vector w0∗ as the initial values for online learning, we have ∇E(w(k))·Δw(k)≤0 if w(k) is updated with the learning rules from (53) to (57).Proof.
For the first case, i.e., if∇E(w(k))·Δw(k-1)<0 holds, then(63)Δwk=-ηk∇Ewk+αkΔwk-1=ηk-∇Ewk+β0ηk-∇Ewk·Δwk-1Δwk-12Δwk-1=ηk-∇Ewk+β0ηk∇EwkΔwk-1cos⁡θkΔwk-12Δwk-1=ηk-∇Ewk+β0ηkcos⁡θk∇EwkΔwk-1Δwk-1,(64)∇Ewk·Δwk=∇ETwkΔwk=ηk-∇Ewk2+β0ηkcos⁡θk∇EwkΔwk-1∇ETwkΔwk-1=-ηk∇Ewk21+β0ηkcos2⁡θk.
As a matter of fact,η(k)≥0 and 1+β0η(k)cos2θ(k)>0 when 0≤β0<1 and 0≤η0≤1, thus ∇E(w(k))·Δw(k)≤0 holds.
For the second case, i.e., if∇E(w(k))·Δw(k-1)≥0, then β(k)=0, α(k)=0, and Δw(k)=-η(k)∇E(w(k)), then(65)∇Ewk·Δwk=∇ETwkΔwk=-ηk∇Ewk2≤0.
The above two cases verify that∇E(w(k))·Δw(k)≤0 holds.
That completes the proof of Theorem3.


## 7. Algorithm for MTN Optimal Control Scheme

Steps of the algorithm for MTN optimal control scheme are summarized as follows.Algorithm 4.
⁡
Step  1. Build the model MTNI for the system (1), train it offline, and take the trained weight parameters as the initial values of MTNI for online identification, where k=0,1,…,N-1. In the process of offline training, M sequence is selected as the persistently exciting input signal for system (1) and model (4).
Step  2. Select an ideal output signal yid(k) for system (1) relative to the given reference signal r(k).
Step  3. Call Algorithm 2 to calculate the optimal control law u∗(k) of the system (1) relative to the ideal output signal yid(k) offline, and the corresponding optimal output signal y∗(k) is taken as the desired output signal.
Step  4. Construct MTNC to fit the optimal control law u∗(k), select a group of initial weight parameter vector w0 in the interval (-1,1) in a random way for offline training, and train MTNC by CG method to obtain the initial weight parameter vector w0∗ for online training.
Step  5. Obtain r(k), u(k), and y(k) by sampling, and take the real-time tracking error e(k) as e(k)=r(k)-y(k), where k=0,1,….
Step  6. Obtain the input signal u(k) of the system (1) by substituting z(k) into MTNC controller (29), and u(k) into (1) and (4) respectively.
Step  7. Utilize the formulas (19)-(23) to adjust the weight parameter vector for MTNI online and the formulas (53)-(57) to tune that for MTNC online.
Step  8. Return to Step 5 to continue the process.


## 8. Simulation Example

To demonstrate the effectiveness of the proposed MTN optimal control scheme, simulation results are presented from the following example, that is, a modification of example 2 derived from [57]. Consider the following SISO nonlinear time-varying discrete system described by the input-output difference equation:(66)yk+1=akyk+bk1+y2k+ckyk-1+uk+cku2k-1,where a(k)=0.96(1-0.5e-0.2k), b(k)=0.1lg(k+1) and c(k)=e-0.2k. The given reference signal is defined as r(k)=1. MTNI takes the 4-15-1 structure with 4 input nodes, 2 powers, and 1 output node, and its input vector is [y(k),y(k-1),u(k-1),u(k)]T. Let z^(k)=[k^1y(k),k^2y(k-1),l^2u(k-1),l^1u(k)]T, and set k^1, k^2, l^1, l^2 as 0.001. As comparison, the system model is also built by the NN identifier (NNI) and PIDNN identifier (PIDNNI). NNI is built from a three-layer NN of the 4-50-1 structure with 4 input nodes, 50 hidden neurons, and 1 output neuron. As confirmed by our simulation, a better identification result can be obtained when 50 hidden neurons are chosen for NNI. The activation functions for the hidden and output layer are chosen as a^h(x)=1/1+e-x and a^o(x)=1-e-x/1+e-x, respectively. PIDNNI is built of a three-layer network of the 4-3-1 structure with 4 input nodes, 3 hidden neurons, and 1 output neuron. The PID neuron structure is obtainable from [58, 59]. At the same time, input vectors for NNI and PIDNNI are the same as for MTNI. For offline identification, the initial weight parameters for MTNI and NNI are chosen randomly in the interval (-1,1). The initial weight parameters between the input and hidden layer for PIDNNI are chosen as w^s1=[1,1,1,1]T, w^s2=[0.1,0.1,-0.1,-0.1]T, and w^s3=[1,1,-1,-1]T. The initial weight parameters between the hidden and output layer for PIDNNI are taken as w^d=[0.1,0.1,0.1]T, and the BP algorithm is chosen for the learning process. Substituting the output y^(k) of MTNI for the output y(k) of the system (66) gives(67)xk=x1k,x2k,x3kT=k^1y^k,k^2y^k-1,l^2uk-1T.And its corresponding extended state space is(68)x1k+1=k^1w^1+k^1w^2x1k+k^1w^3x2k+k^1w^4x3k+k^1l^1w^5uk+k^1w^6x12k+k^1w^7x1kx2k+k^1w^8x1kx3k+k^1l^1w^9x1kuk+k^1w^10x22k+k^1w^11x2kx3k+k^1l^1w^12x2kuk+k^1w^13x32k+k^1l^1w^14x3kuk+k^1l^12w^15u2k,x2k+1=k^2k^1x1k,x3k+1=l^2uk,(69)y^k=1k^1x1kwhere w^j (j=1,2,…,15) can be obtained offline.To generate the control input, MTNC takes the 3-10-1 structure with 3 input nodes, 2 powers, and 1 output node, its input vector beingz(k)=[k1u(k-1),l1e(k-1),l2e(k-2)]T, where k1,l1,l2 are set as 0.001. As comparison, the NN adaptive controller (NNAC) is built from a three-layer NN of the 3-50-1 structure with 3 input neurons, 50 hidden neurons, and 1 output neuron. As demonstrated by our experiments, a better control result can be obtained when 50 hidden neurons are chosen for NNAC. The activation functions for the hidden and output layer are chosen as ah(x)=1/1+e-x and ao(x)=1-e-x/1+e-x, respectively. The PIDNN adaptive controller (PIDNNAC) takes the 2-3-1 structure with 2 input nodes, 3 hidden neurons, and 1 output node, respectively. Before the real-time control process, the MTNC controller is generated automatically by fitting the optimal control law u∗(k) using the CG method, with the initial weight values of MTNC chosen randomly in the interval (-1,1) and the iterations of fitting set as ite=100. For online identification and control, the initial weight values of MTNI and MTNC are taken as those from the offline learning. The novel adaptive BP algorithm with the learning rules from (19) to (23) is employed to update the weight values of MTNI to build the system model in real time, and the formulas from (53) to (57) are utilized to update the those of MTNC to implement the output tracking control of the system (66) relative to the given reference signal r(k). The initial weight values for NNAC are chosen randomly in the interval (-1,1), those of PIDNNAC between the input and hidden layer are set as ws1=[1,1]T, ws2=[0.1,-0.1]T, and ws3=[1,-0.1]T, and those of PIDNNI are set as wd=[0.1,0.1,0.1]T. BP algorithm is adopted to train NNAC and PIDNNAC online. For online identification and control, the initial values of η^(0), η^0, β^0, and η^max for MTNI are set as 0.2, 0.001, 0.01, and 0.5, respectively, and those of η(0), η0, β0, and ηmax for MTNC as 0.2, 0.001, 0.01, and 0.5, respectively. The learning factor is set as α^N = 0.2 for NNI, αN = 0.2 for NAC, α^P=0.3 for PIDNNI, and αP=0.3 for PIDNNAC. The tracking results and errors for the system with MTNC, NNAC, and PIDNNAC are shown as in Figures 9 and 10, respectively, and the corresponding control inputs are presented by Figure 11.Figure 9
The tracking trajectories.Figure 10
The tracking errors.Figure 11
The control inputs.Remark 5.
In Figures9, 10, and 11, r represents the given reference signal; yMC, yNAC, and yPNAC are the actual output responses for MTNC, NN, and PIDNN control schemes, respectively; eMC, eNAC, and ePNAC represent the corresponding tracking errors; uMC, uNAC, and uPNAC are the corresponding control inputs.
From Figures9 and 10, it can be seen that the overshoots are 31.19%, 88.43%, and 67.29%; the performance index isE<10-3 after the iterations of 14, 20, and 24 with MTN, NN, and PIDNN control schemes; the steady-state error worked out by taking the average of the absolute errors from the iterations of 14 to 300 is 0.0024 with the MTN optimal control scheme, 0.0011 taken from 20 to 300 with the traditional NN adaptive control scheme, and 8.0383×10-4 taken from 24 to 300 with the PIDNN control scheme. Simulation results demonstrate that the proposed control scheme outperforms the others.
(a) Noise Interference Experiments. At time 100, a Gaussian white noise with the mean of 0 and the standard deviation of 0.2 is added to system (66). The simulation results are presented in Figures 12 and 13, and the corresponding control inputs are shown in Figure 14.
As demonstrated by Figures12 and 13, the MTN optimal control scheme promises a better robustness for the noise interference than the other two.
For clearer illustration of the robustness of the proposed control scheme, the noise interference is expanded 30 times based on the above discussion, and the simulation results are presented in Figures15 and 16 and the corresponding control inputs are shown in Figure 17.
(b) Input Superposition Experiments. When the external superposition d(k)=0.05sin(0.005πk) is added to the given reference signal r(k), the simulation results are shown in Figures 18 and 19 and the corresponding control inputs are shown as in Figure 20.Figure 12
The tracking trajectories with noise interference.Figure 13
The tracking errors with noise interference.Figure 14
The control inputs with noise interference.Figure 15
The tracking trajectories with expanded noise interference.Figure 16
The tracking errors with expanded noise interference.Figure 17
The control inputs with expanded noise interference.Figure 18
The tracking trajectories with signal superposition.Figure 19
The tracking errors with signal superposition.Figure 20
The control inputs with signal superposition.As demonstrated by Figures9–20, the proposed MTN scheme is valid and of desirable for the alteration of the system parameters.Remark 6.
For each control period, 29 times of multiplication operations and 14 times of addition operations are needed for MTNI, and 19 times of multiplication and 9 times of addition for MTNC. That is, 48 times of multiplication and 23 times of addition are required by the MTN scheme. However, with the exponential function expanded into finite terms with 2 powers, 8-time multiplication and 9-time addition operations are needed for each hidden node, and 57-time multiplication and 61-time addition operations are required for the output layer node of NNI. Meanwhile, data normalization and inverse normalization are considered for the online identification process, and the two functions are defined asfo(x)=2/πarctanx and fi(x)=tanπ/2x, respectively. For each normalization, 5-time multiplication and 1-time addition operations are needed with fo(x) expanded into finite terms with 3 powers; 5-time multiplication and 1-time addition operations are required with fi(x) expanded into finite terms with 3 powers. For NNAC, 1-time multiplication and 3-time addition operations are demanded for the input layer, 7-time multiplication and 8-time addition operations for each hidden node, and 57-time multiplication and 61-time addition operations for each output layer node. That is, for each control period, 875-time multiplication and 977-time addition operations are needed implementing the traditional neural network adaptive control scheme. With the PIDNN control scheme, for PIDNNI, 12-time multiplication operations and 11-time addition operations are required for the hidden layer nodes, and 3-time multiplication and 2-time addition operations for the output layer node. Meanwhile, data normalization and inverse normalization are also considered for the online identification process, and the two functions are defined as fo(x)=2/πarctanx and fi(x)=tanπ/2x, respectively. Then, for each normalization, 5-time multiplication and 1-time addition operations are needed with fo(x) expanded into finite terms with 3 powers; 5-time multiplication and 1-time addition operations are required with fi(x) expanded into finite terms with 3 powers. For PIDNNAC, 6-time multiplication and 5-time addition operations are demanded for the hidden layer nodes, and 3-time multiplication and 2-time addition operations for the output layer node. That is, for each control period, 34-time multiplication and 22-time addition operations are needed implementing the PIDNN control scheme, while fewer operations are required by the MTN and PIDNN schemes. And compared to the type TMS320F28335 of DSP with the dominant frequency of 150MHz, the computation time for each control period is 473.3ns, 12346.7ns, and 373.3ns with the MTN, NN, and PIDNN control scheme, respectively. As verified by the simulation and calculation results, the MTN and PIDNN control schemes promise more desirable real-time control than the NN scheme.


## 9. Conclusions

For the SISO nonlinear time-varying discrete system without mechanism model, a MTN optimal control scheme has been proposed to secure its real-time output tracking control based on the given reference signal.Main contributions of the paper can be summarized as follows:(1) MTN optimal control scheme has been proposed for general nonlinear time-varying discrete system control design; (2) MTNI identifier and MTNC controller have been built to simplify the network structure and raise the convergence speed; (3) initial value selection scheme for the weight parameters of the controller has been developed; (4) novel adaptive BP learning algorithm has been proposed to adjust the weight parameters for a faster convergence speed.Simulation results show that the proposed control scheme is effective and capable of enabling the system’s actual output response to well track the given reference signal in real time.



---

*Source: 3219140-2018-07-03.xml*