Dynamical system identification using radial-basis

function network

Stanislaw H. Żak ^{∗}

February 13, 1999

Abstract

In this section, we present a radial-basis function (RBF) neural network identifier of an unknown nonlinear dynamical system using only the system’s input and output data. Our exposition is based on the paper by Liu, Kadirkamanathan, and Billings [?]. The neural identifier is an adaptive system for which the learning algorithm is developed using Lyapunov techniques. The algorithm allocates resources, by adding new nodes, to obtain the desired identification accuracy.

1 Radial-Basis Function Networks

A diagram of a radial-basis function (RBF) network is shown in Figure 1. The functions ϕ (∥x − x^{(i)}∥), i = 1, 2, …, N are called the radial-basis functions, where ∥ · ∥ denotes a p-norm that is usually the Euclidean 2-norm. The real paramters w_{i}, i = 1, 2, …, N are the weights. A common choice for ϕ is

ϕ (∥x − x^{(i)}∥) = exp (−\frac{∥x − x^{(i)}∥}{2σ_{i}^{2}} ) ,  (1)

∗School of Electrical and Computer Engineering, Purdue University, West Lafayette, IN 47907

Figure 1: A radial-basis function (RGF) network.

where x^{(i)} is the center of the basis function and σ_{i} is sometimes called the radius of the radial basis function. The function given by (1) is called the Gaussian radial basis function.

The RBF networks are used to perform interpolation. The interpolation problem can be formulated as: For a given set of N different points,

{x^{(i)} : x^{(i)} ∈ ℝ^{n}, i = 1, 2, …, N}

and a corresponding set of N real numbers

{d_{i} : d_{i} ∈ ℝ, i = 1, 2, …, N}

find a function

F : ℝ^{n} → ℝ

such that

F (x^{(i)}) = d_{i}, i = 1, 2, …, N.  (2)

The function that the RBF network implements is

F (x) = ∑_{i=1}^{N} w_{i}ϕ (∥x − x^{(i)}∥) .  (3)

\begin{matrix}x^{(i)} & d_{i} \\ 0 & 0.8 \\ 6 & 0.2 \\ 7 & 0.4 \\ 8 & 0.6 \\ 10 & 0.3 \end{matrix}

Substituting the interpolation conditions given by (2) into (3) gives

[ \begin{matrix}ϕ_{11} & ϕ_{12} & ⋯ & ϕ_{1N} \\ ϕ_{21} & ϕ_{22} & ⋯ & ϕ_{2N} \\ ⋮ & ⋮ & ... & ⋮ \\ ϕ_{11} & ϕ_{12} & ⋯ & ϕ_{1N} \end{matrix} ] [ w_{1} ; w_{2} ; ⋮ ; w_{N} ] = [ d_{1} ; d_{2} ; ⋮ ; d_{N} ] .  (4)

Let w = [ w_{1} w_{2} ⋯ w_{N} ], d = [ d_{1} d_{2} ⋯ d_{N} ], and

Φ = [ \begin{matrix}ϕ_{11} & ϕ_{12} & ⋯ & ϕ_{1N} \\ ϕ_{21} & ϕ_{22} & ⋯ & ϕ_{2N} \\ ⋮ & ⋮ & ... & ⋮ \\ ϕ_{11} & ϕ_{12} & ⋯ & ϕ_{1N} \end{matrix} ] .

Then, using the above notation, we represent (4) as

Φw = d.  (5)

The matrix Φ is called the interpolation matrix. If the interpolation matrix is nonsingular, then we can solve (5) for w to obtain

w = Φ^{−1}d.  (6)

Example 1 We construct an RBF network for the interpolation data given in Table 1. We set σ_{i} = 0, i = 1 …, 5 and evaluate (6) to obtain

w = [ 0.8000 0.1459 −0.0390 0.5737 0.2227 ]^{T}

Figure 2: Results of interpolating using radial-basis functions.

The interpolation result is depicted in Figure 2. In Figure (2)(a), we used solid line to show the result of interpolating using the Gaussian radial-radial basis functions, while broken line indicates the result of straight line interpolation. The plots of the basis functions that contribute to the interpolation are shown in Figure (2)(b).

2 Identification of a single-input single-state dynami-

cal system

We first discuss the identification of the single-input single-state dynamical system using neural network. Let the dynamical system be modeled by

ẋ(t) = f(x(t), u(t)), x(0) = x_{0},  (7)

where f : ℝ × ℝ → ℝ is an unknown, possibly nonlinear, function of x and u to be identified by a neural network. Adding and subtracting ax(t), where a > 0 is a design parameter, to the right-hand side of (7) gives

ẋ(t) = −ax(t) + ax(t) + f(x(t), u(t))
= −ax(t) + g(x(t), u(t)),  (8)

where g(x(t), u(t)) = ax(t) + f(x(t), u(t)). The reason for the above manipulations will become apparent when we discuss the learning rules. We will now proceed with the construction of a neural identifier of the dynamical system (7). The identifier will be described by

˙x̂(t) = −ax̂(t) + ĝ(x(t), u(t), m_{k}, r_{k}), x̂(0) = x̂_{0},  (9)

where x̂ is the state of the identification network, m_{k} and r_{k} are the adjustable parameters of the identifier. The function ĝ has the form,

ĝ(x(t), u(t), m_{k}, r_{k}) = ∑_{k=1}^{K} w_{k} exp (−\frac{1}{r_{k}^{2}} ((x(t) − m_{xk})^{2} + (u(t) − m_{uk})^{2})) ,  (10)

where m_{k} = [ m_{xk} m_{uk} ], k = 1, 2, …, K. For the sake of brevity, we denote

ϕ_{k}(x, u, m_{k}, r_{k}) = exp (−\frac{1}{r_{k}^{2}} ((x(t) − m_{xk})^{2} + (u(t) − m_{uk})^{2})) .

Using the above notation, we represent (10) as

ĝ(x(t), u(t), m_{k}, r_{k}) = ∑_{k=1}^{K} w_{k}ϕ_{k}(x(t), u(t), m_{k}, r_{k})  (11)

Suppose now that the basis functions ϕ_{k}s are given to us. The problem of choosing the basis functions will be discussed later. Suppose also that we know the “optimal” weights w_{k}^{∗}, k = 1, 2, …, K, that result in the least identification error between g in (8) and ĝ in (9), where the identification error is defined as

e(t) = g(x(t), u(t)) − ∑_{k=1}^{K} w_{k}^{∗}ϕ_{k}(x(t), u(t), m_{k}, r_{k})  (12)

Figure 3: A configuration for the radial-basis function neural identifier.

Using the above notation, we represent the dynamical system model (8) as

ẋ(t) = −ax(t) + ∑_{k=1}^{K} w_{k}^{∗}ϕ_{k}(x(t), u(t), m_{k}, r_{k}) + e(t).  (13)

The neural identifier, on the other hand can be represented as

˙x̂(t) = −ax̂(t) + ∑_{k=1}^{K} w_{k}ϕ_{k}(x(t), u(t), m_{k}, r_{k})  (14)

A schematic of the above identification architecture is depicted in Figure 3.

3 Learning algorithm for the RBF neural identifier

We assume that both the state x and the input u of the dynamical system model (7), and hence that of (8), are bounded. We next define the state estimation error to be

e_{x} = x − x̂,

while the weights’ estimation errors, ξ_{k} = w_{k}^{∗} − w_{k}, k = 1, …, K.

We assmue that the state estimation error, e_{x}, is bounded by some finite constant e_{K}, that is, |e_{x}| ≤ e_{K}. However, e_{K} is unkown to us. Next, we use (13) and (14) to form the differential equation describing the state error dynamics. We have

ė_{x} = ẋ − ˙ x̂
= −ax + ∑_{k=1}^{K} w_{k}^{∗}ϕ_{k}(x, u, m_{k}, r_{k}) + e + ax̂ − ∑_{k=1}^{K} w_{k}ϕ_{k}(x, u, m_{k}, r_{k})
= −ae_{x} + ∑_{k=1}^{K} ξ_{k}ϕ_{k}(x, u, m_{k}, r_{k}) + e.  (15)

Our goal now is to use the second method of Lyapunov to devise a learning law for the RBF identifier that would yield stable identification scheme in the sense that the state error e_{x} decreases as the time increases and the weights do not diverge. For this purpose, we will use the following Lyapunov function candidate,

V (e_{x}, ξ) = \frac{1}{2} (e_{x}^{2} + \frac{1}{α} ∑_{k=1}^{K} ξ_{k}^{2}) ,  (16)

where α > 0 is a design parameter and ξ = [ ξ_{1} ξ_{2} ⋯ ξ_{K} ]^{T}. We now evaluate the time derivative of V on the trajectories of (15),

V̇ = e_{x}ė_{x} + \frac{1}{α} ∑_{k=1}^{K} ξ_{k}ξ̇_{k}
= e_{x} (−ae_{x} + ∑_{k=1}^{K} ξ_{k}ϕ_{k}(x, u, m_{k}, r_{k}) + e) + \frac{1}{α} ∑_{k=1}^{K} ξ_{k}ξ̇_{k}
= −ae_{x}^{2} + \frac{1}{α} ∑_{k=1}^{K} ξ_{k} (αe_{x}ϕ_{k}(x, u, m_{k}, r_{k}) + ξ̇_{k}) + e_{x}e.  (17)

Let

ξ̇_{k} = −αe_{x}ϕ_{k}(x, u, m_{k}, r_{k}).  (18)

Then (17) simplifies to

V̇ = −ae_{x}^{2} + e_{x}e.  (19)

Since ẇ_{k} = −ξ̇_{k}, we can represent (18) as

ẇ_{k} = αe_{x}ϕ_{k}(x, u, m_{k}, r_{k})  (20)

We further evaluate (19) to obtain

V̇ ≤ −ae_{x}^{2} + |e_{x}| e_{K}
= −a |e_{x}| (|e_{x}| − \frac{e_{K}}{a} )  (21)

It follows from the above that if there is no state estimation error, that is, e_{K} = 0, then V̇ is negative semi-definite in the [ e_{x} ξ^{T} ]^{T} -space. In such a case the origin of the [ e_{x} ξ^{T} ]^{T} - space is stable, which guarantees boundedness of the weights w_{k}s. In the presence of the state estimation error, if |e_{x}| < e_{K}/a, then it is possible that V̇ > 0, which implies that the weights w_{k}s may drift to infinity with time. To avoid this problem, Liu, Kadirkamanathan, and Billings [?] propose to set ẇ_{k} = 0 if |e_{x}| < e_{K}/a. Also since e_{K} is unknown, a bound is imposed on the weights as well as on the error. The learning rule becomes

ẇ_{k} = \{ \begin{matrix}αe_{x}ϕ_{k}(x, u, m_{k}, r_{k}) & if & |e_{x}| ≥ e_{0} and ∥w∥ ≤ M \\ 0 & & otherwise \end{matrix} (22)

where the bounds e_{0} and M are design parameters, and w = [ w_{1} ⋯ w_{K} ]^{T}. With the above modified learning law, whenever |e_{x}| ≥ e_{0} ≥ e_{K}/a, the Lyapunov derivative V̇ is always semi-definite. However, when e_{0} ≤ |e_{x}| ≤ e_{K}/a, the weights may increase with time because it is possible that V̇ > 0. But the weights, as well as the estimation error, e_{x}, are upper-bounded. When any of the two bounds is reached, more basis functions maybe needed to improve the approximation of g. New basis functions are added to reduce the error until V̇ < 0. We now describe an algorithm for adding radial-basis functions to the identification network.

4 Growing radial-basis function network

As we mentioned above, when V̇ > 0, the identification scheme is no longer stable, but the learning rule (22) guarantees that the adjustable parameters are bounded. Once the parameter adjustments can no longer reduce the error and the bounds are reached, additional basis functions in the RBF network should be added. This continues until the approximation leads to the error reduction so that V̇ < 0.

In the proposed identification scheme, we use the gaussian radial-basis functions, whose areas of “influence” are determined by the functions’ centers and widths. As the function to be approximated moves further away from a particular basis function, the contribution of that basis function decreases. This is due to the exponential form of the gaussian basis functions. As z increases exp(−z) decreases. The maximum contribution of a particular basis function to approximating a given function occurs when the center of the basis function and the approximated function are at the same point, since in this case we have exp(0) = 1. This heuristic argument leads to the following rule for determining the centers of new basis functions to be added at a time instance t,

m_{K+1} = [ round (x(t)/δ_{x}) δ_{x} round (u(t)/δ_{u}) δ_{u} ]^{T} ,  (23)

where the round operator rounds the elements on which it operates to the nearest integer. For example, let x = [ −1.8 −0.3 3.7 ]. Then, round(x) = [ −2 0 4 ]. The parameters, δ_{x} and δ_{u} represent the spacing of the centers of the basis functions. We can think as having a grid over the input-state space where the spacing between the vertices of the grid are determined by the parameters δ_{x} and δ_{u}. Each vertex of the grid represents the potential center of a basis function. The radius of the new basis function and the adjusted radii are determined from the relation,

r_{k} = κ ∥m_{k} − m_{k}^{∗}∥ ,  (24)

where κ > 0 is a design parameter and

m_{k}^{∗} = arg min_{{i=1,...,K \atop m_{i}≠m_{k}}} \{∥m_{k} − m_{k}^{∗}∥\} .  (25)

Thus, m_{k}^{∗} is the nearest center of all basis functions to the center of the k-th basis function. With the above background, we can now present the rules for adding a new basis function to the RBF network. Initially, the network contains no basis functions. Then, once the following two conditions are satisfied, a new basis function is added. The first condition is:

min_{k=1,...,K} |x(t) − m_{xk}| > δ_{x}/2 or min_{k=1,...,K} |x(t) − m_{uk}| > δ_{u}/2.  (26)

The above will be satisfied if the state or input, at time t, is closer to an available grid vertex than to any of the existing centers of the basis functions. The basis function at the new vertex would be the closest to the existing state or input and would have the greatest influence in the approximation of g. The second condition is

|x(t) − x̂(t)| > e_{max},  (27)

where e_{max} is the maximum allowable state estimation error. If condition (27) is satisfied, a new basis function is added to reduce the state estimation error.

Example 2 The above technique for identification of dynamical systems using growing radial-basis function network is illustrated with the following system,

ẋ = 1.1 (1 − x + 2xu − x^{2}) exp (−0.5x^{2} − u^{2}) , x(0) = 0,  (28)

where u = cos(t/100). In our simulations, we discretized the above differential equation using the Euler method with the step size dt = 0.01 sec. The initial weights and the initial state of the identifier were set to zero. We use the following parameter values: e_{0} = 0.001, e_{max} = 0.005, α = 5, κ = 3. The MATLAB code appears below.

t=1;dt=0.01;w(1)=0;
i=0;g(1)=0;xhat(1)=0;
emax=0.005;e0=0.001;
dx=0.05;du=0.05;temp=0;
alpha=5;kappa=3;
r(1)=kappa*dx;
x(1)=0;
m1(1)=x(1);
m2(1)=cos(0);
tfinal=1000;
while t<tfinal

u(t)=cos(t*dt);
xdot(t)=1.1*(1-x(t)+2*x(t)*u(t)-x(t)^2)*exp(-0.5*x(t)^2-u(t)^2);
x(t+1)=x(t)+xdot(t)*dt;
for z=1:100
g(t)=0;
if i>0
   for n=1:i
     g(t)=g(t)+w(n)*exp(-((x(t)-m1(n))^2+(u(t)-m2(n))^2)/(r(n)^2));
   end %end for
end %end if
xhatdot(t)=-0.5*xhat(t)+g(t);
xhat(t+1)=xhat(t)+xhatdot(t)*dt;
e(t)=x(t)-xhat(t);
xmin=abs(x(t)-m1(1));
umin=abs(u(t)-m2(1));
if i>0
   if e0 <= abs(e(t))
       if norm(w) <= 1.5*sqrt(i)
           for n=1:i
    wdot=alpha*e(t)*exp(-((x(t)-m1(n))^2+(u(t)-m2(n))^2)/(r(n)^2));
              w(n)=w(n)+wdot*dt;
              if abs(x(t)-m1(n))<xmin
                   xmin=abs(x(t)-m1(n));
              end %end if
              if abs(u(t)-m2(n))<umin
                   umin=abs(u(t)-m2(n));
              end
           end

    end
end
end
if abs(e(t))>emax
if xmin > (0.5*dx)
    temp=i+1;
end
if umin > (0.5*du)
    temp=i+1;
end
end
if i<temp
i=temp;
m1(i)=round(x(t)/dx)*dx;
m2(i)=round(u(t)/du)*du;
w(i)=0;
r(i)=3*dx;
for a=1:i
    m1min=0;
    m2min=0;
      if a~=i
           mk=sqrt((m1(a)-m1(a+1))^2+(m2(a)-m2(a+1))^2);
           m1min=m1(a+1);
           m2min=m2(a+1);
      elseif i~=1
           mk=sqrt((m1(a)-m1(a-1))^2+(m2(a)-m2(a-1))^2);
           m1min=m1(a-1);
           m2min=m2(a-1);

                 end
     for b=1:i
       if a~=b
         if mk>sqrt((m1(a)-m1(b))^2+(m2(a)-m2(b))^2)
             mk=sqrt((m1(a)-m1(b))^2+(m2(a)-m2(b))^2);
             m1min=m1(b);
             m2min=m2(b);
         end
       end
      end
       r(a)=kappa*sqrt((m1(a)-m1min)^2+(m2(a)-m2min)^2);
      end
    w(i)=dt*alpha*e(t)*exp(-((x(t)-m1(i))^2+(u(t)-m2(i))^2)/(r(i)^2));
   end
end
t=t+1;
end
plot(1:tfinal,x,1:tfinal,xhat)
xlabel('Iteration','Fontsize',16)
ylabel('x and its estimate','Fontsize',16)

We mention that the z learning loop from 1 to 100 was added to allow the approximation to learn out the error. Using the above program, we generated a number of plots illustrating the performance of the RBF neural identifier. The last three lines in the MATLAB code show how we obtained the plots of Figure 4. The plots in the remaining figures were obtained using similar commands. In Figure 4, plots of the true state and its estimate versus time are shown. A plot of the error, e_{x}, between the true state and its estimate is shown in Figure 5. Initially, the error exceeds the bound e_{max}. Adding additional basis functions drives the error back down. The plots shown in Figure 6 depict the true state and its estimate versus the Figure 4: Plots of the true state and its estimate versus time. The estimated state “hugs” closely the actual state.

Figure 5: A plot of the error, e_{x}, between the true state and its estimate.

Figure 6: A plot of the true and its estimate versus the input u.

input u. Finally, in Figure 7, plots of the derivative of the true state versus the true state and the derivative of the estimated state versus the state estimate are shown.

5 Identification of multivariable systems

We now extend the above concepts to multi-input multi-state dynamical systems. In our derivation of the update law for the weights of the RBF network, we use a result from matrix algebra that we present next. Recall the definition of the trace of a square matrix as given by (??). If Γ ∈ ℝ^{m×n}, then by (??) its Frobenius norm satisfies

∥Γ∥_{F}^{2} = trace (Γ^{T} Γ) = ∑_{i=1}^{m} ∑_{j=1}^{n} γ_{ij}^{2}.

If the elements, γ_{ij}, of Γ are functions of time, then

\frac{d}{dt} (trace (Γ^{T} Γ)) = 2 ∑_{i=1}^{m} ∑_{j=1}^{n} γ_{ij}γ̇_{ij} = 2 trace (Γ^{T} Γ̇) ,  (29)

Figure 7: Plots of the derivative of the true state versus the true state and the derivative of the estimated state versus the state estimate.

where Γ̇ = \frac{dΓ}{dt}. With the above background, we can proceed with the anlaysis of the RBF neural identifier for the multivariable systems.

The multivariable dynamical systems that we consider are modeled by

ẋ(t) = f(x(t), u(t)), x(0) = x_{0},  (30)

where u(t) ∈ ℝ^{m}, x(t) ∈ ℝ^{n}, and f : ℝ^{n} × ℝ^{m} → ℝ^{n} is the nonlinear function to be identified. Let A ∈ ℝ^{n×n} be a Hurwitz matrix, that is, the matrix whose eigenvalues are all in the right-half palne. Since A, then by the theorem of Lyapunov, which is Theorem ??, for any given real symmetric positive definite (r.s.p.d.) matrix Q the solution P to the continuous Lyapunov matrix equation,

A^{T} P + P A = −Q  (31)

is also (real symmetric) positive definite. We use this fact later when we deal with a learning law for the RBF neural identifier.

We represent the model (30) in an equivalent form,

ẋ(t) = Ax(t) + g(x(t), u(t)), x(0) = x_{0},  (32)

where

g(x(t), u(t)) = −Ax(t) + f(x(t), u(t)).

Let the k-th basis function have the form

ϕ_{k}(x, u) = exp (−\frac{1}{r_{k}^{2}} ‖[ \begin{matrix}x \\ u \end{matrix} ] − m_{k}‖^{2}) ,

and let

Φ_{K} = [ ϕ_{1}(x, u) ϕ_{2}(x, u) ⋯ ϕ_{K}(x, u) ]^{T}

and

ĝ(x, u, m_{k}, r_{k}) = W_{K}Φ_{K},

where W_{K} is the n × K weight matrix and m_{k} = [ m_{xk}^{T} m_{uk}^{T} ]^{T} is the center of the k-th radial basis function in the [ x^{T} u^{T} ]^{T} -space. The RBF neural identifier has the form,

˙x̂(t) = Ax̂(t) + ĝ(x, u, m_{k}, r_{k}), x̂(0) = x_{0}.  (33)

Let W _{K}^{∗} be the optimal weight matrix that minimizes the modeling error,

e(t) = g(x(t), u(t)) − W _{K}^{∗}Φ_{K}(x(t), u(t)).

Using the above definition, we represent the system model (32) as

ẋ(t) = Ax(t) + W _{K}^{∗}Φ_{K}(x(t), u(t)) + e(t).  (34)

We define the state estimation error and the weight estimation error as ex = x − x̂ and Γ_{K} = W _{K}^{∗} − W_{K}.

Next, we form the differential equation describing the state estimation error dynamics. We have

ėx = ẋ − ˙ x̂
= Ax + W _{K}^{∗}Φ_{K}(x, u) + e − Ax̂ − ĝ(x, u, m_{k}, r_{k})
= A (x − x̂) + W _{K}^{∗}Φ_{K}(x, u) − W_{K}Φ_{K}(x, u) + e
= Aex + Γ_{K}Φ_{K}(x, u) + e.  (35)

We use the following Lyapunov function candidate,

V = V (ex, Γ_{K}) = \frac{1}{2} (ex^{T}P ex + \frac{1}{α}trace (Γ_{K}Γ_{K}^{T})) .  (36)

Using (29), we evaluate the time derivative of V on the solutions of the state estimation error equation to get

V̇ = ex^{T}P ėx + \frac{1}{α}trace (Γ_{K} Γ̇_{K}^{T})
= ex^{T} (P Aex + P Γ_{K}Φ_{K}(x, u) + P e) + \frac{1}{α}trace (Γ_{K} Γ̇_{K}^{T})
= \frac{1}{2}ex^{T} (A^{T} P + P A) ex + ex^{T}P Γ_{K}Φ_{K}(x, u) + \frac{1}{α}trace (Γ_{K} Γ̇_{K}^{T}) + ex^{T}P e. (37)

Note that by (??) for non-square matrices, we have

ex^{T}P Γ_{K}Φ_{K}(x, u) = trace (ex^{T}P Γ_{K}Φ_{K}(x, u))
= trace (Γ_{K}Φ_{K}(x, u)ex^{T}P ) .  (38)

Observe that

Ẇ_{K} = −Γ̇_{K}

Using the above and (38) along with the Lyapunov equation (31), we write

V̇ = −ex^{T}Qex + \frac{1}{α} (trace (αΓ_{K}Φ_{K}(x, u)ex^{T}P ) − trace (Γ_{K} Ẇ _{K}^{T})) + ex^{T}P e.  (39)

We propose to choose the following learning rule,

Ẇ_{K} = αP exΦ_{K}^{T}(x, u)  (40)

With the above weight matrix update law, the Lyapunov derivative (39) evaluates to

V̇ = −ex^{T}Qex + ex^{T}P e
≤ −λ_{min}(Q)∥ex∥^{2} + λ_{max}(P)∥e∥∥ex∥
= −λ_{min}(Q) (∥ex∥ − \frac{λ_{max}(P)}{λ_{min}(Q)}∥e∥) ∥ex∥.  (41)

It follows from the above that if the magnitude of the state estimation error is such that ∥ex∥ < \frac{λ_{max}(P)}{λ_{min}(Q)}∥e∥, then it may happen that V̇ > 0. This in turn may lead to the weights’ drift to infinity over time. To prevent this, we modify (40) following the analysis of the single-input single-state case. The modified weight update law is

Ẇ_{K} = \{ \begin{matrix}αP exΦ_{K}^{T}(x, u) & if & ∥ex∥ ≥ e_{0} and ∥W_{K}∥ ≤ M \\ O & & otherwise \end{matrix} (42)