# Twin Support Vector Machine Method for Identification of Wiener Models

Authors: Mujahed Al-Dhaifallah
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Twin support vector regression is applied to identify nonlinear
Wiener system, consisting of a linear dynamic block in series with static nonlinearity.
The linear block is expanded in terms of basis functions, such as Laguerre or Kautz
filters, and the static nonlinear block is determined using twin support vector machine
regression. Simulation of a control valve model and pH neutralization process have been
presented to show the features of the proposed algorithm over support vector machine
based algorithm.


---

## Body

## 1. Introduction

The goal of system identification is to find a model, within a selected class of models, that produces the best predictions of a system’s output. In general, one forms a cost function that depends on some norm of the prediction errors and finds the model that minimizes this cost function. Since the model is an approximation to the true system, there is a trade-off between the complexity of the structure of the model and the accuracy of its predictions. In many cases, linear models can be used to produce accurate predictions of a system behavior, particularly, if it is restricted to operating within a narrow region. However, if the model is required to cover a broader operating region, then a nonlinear model may be required [1]. Block structured models, cascades of static nonlinearities and dynamic linear systems, are often a good trade-off as they can represent some dynamic nonlinear systems very accurately but are nonetheless quite simple. Common nonlinear models of this type are the Wiener and Hammerstein models [2]. Many algorithms have been proposed to identify Wiener models [3–11]. As one might notice from these researches, the extensive knowledge about linear time invariant (LTI) system representations was applied to the dynamic linear blocks. On the other hand, finding an effective representation for the nonlinearity is an active area of research. Traditionally, the nonlinearity is represented by a polynomial because it is simple and easy to estimate. However, the problem with polynomial approximation is that it cannot deal with many common nonlinearities (saturation, threshold, dead zone, etc.). Better approximation of these nonlinearities can be achieved by spline functions. However, spline functions are defined by a series of knot points which must either be chosen a priori or treated as model parameters and included in the (nonconvex) optimization. Neural networks are another tool to approximate nonlinear functions. Their powerful approximation abilities make them attractive. However, the need to specify the neural network topology in terms of the number of nodes and layers and the need to solve nonconvex optimization complicate their implementation. Recently, support vector machines (SVMs) and least squares support vector machines (LS-SVMs) have demonstrated powerful abilities in approximating linear and nonlinear functions [12, 13]. In contrast with other approximation methods, SVMs do not require a priori structural information. Furthermore, there are well established methods with guaranteed convergence (ordinary least squares, quadratic programming) for fitting LS-SVMs and SVMs [14]. One of the common drawbacks of SVM based algorithms is that they are computationally heavy. Recently, a twin support vector machine (TSVM) regression algorithm has been proposed [15]. The formulation of TSVM is very close to classical SVM except that it aims to enclose data points between two parallel planes such that each plane is closer to one class and is as far as possible from the other. Such formulation reduces the computation complexity which makes the TSVM one of the common methods in machine learning. Lately, many extensions of TSVM classifier have been proposed to improve its performance in certain aspects. A TSVM classifier has been extended to put the TSVM to multitask learning [16]. A generalized framework of TSVM for learning from labeled and unlabeled data was investigated by [17, 18].Recently, Tötterman and Toivonen [19] have developed a new algorithm to identify Wiener models based on SVM regression, where the linear part is described by a basis filter expansion while the nonlinear part is represented by SVM. In this work, TSVM regression is used to formulate an identification algorithm for Wiener models. Simulation examples are presented to show the virtues of the proposed algorithm over Tötterman’s algorithm. The outline of this paper is as follows: TSVM theory is reviewed in Section 2. In Section 3, an algorithm for the identification of Wiener models based on twin support vector machine is proposed. Section 4 presents two illustrative examples to test the proposed algorithm. In Section 5, concluding remarks are given.


## 2. Twin Support Vector Machines Regression

Twin support vector regression (TSVR) is obtained by solving the following pair of quadratic programming problems (QPPs):(1)min12Y-eε1-KA,ATw1-eb1T×Y-eε1-KA,ATw1-eb1+C1ε1eTξs.t.Y-KA,ATw1-eb1⩾eε1-ξ,ξ⩾0,(2)min12Y-eε2-KA,ATw2-eb2T×Y-eε2-KA,ATw2-eb2+C2ε2eTηs.t.KA,ATw2+eb2-Y⩾eε2-η,η⩾0,where C1,C2>0, ε1,ε2⩾0 are parameters, ξ, η are slack variables, e is vector of ones, and K·,· is a nonlinear kernel function. Given training data points A,Y where A and Y represent input and output vectors, respectively, the TSVR algorithm finds two functions ϕ1x=w1TKA,x+b1, which determines the ε1-insensitive down bound regressor, and ϕ2x=w2TKA,x+b2, which determines the ε2-insensitive up bound regressor. The end regressor is computed as the mean of these two functions. The geometric interpretation is given in Figure 1. The objective function of (1) or (2) is the sum of the squared distances from the shifted functions y=w1TKA,x+b1+ε1 or y=w2TKA,x+b2-ε2 to the training points. Therefore, minimizing it leads to the function ϕ1x or ϕ2x. The constraints require the estimated function ϕ1x or ϕ2x to be at a distance of at least ε1 or ε2 from the training points. That is, the training points should be larger than the function ϕ1x at least ε1, while they should be smaller than the function ϕ2x at least ε2. The slack variables ξ and η are introduced to measure the error whenever the distance is closer than ε1 or ε2. The second term of the objective function minimizes the sum error variables, thus attempting to overfit the training points.Figure 1
The geometric interpretation of TSVR.To find the solution of (1) and (2), the dual problems are needed to be derived. The optimization problems (1) and (2) just described are the primal problems for regression. To formulate the corresponding dual problems, the Lagrangian function L is written. Then, L is minimized with respect to the weight w and slack variables ξ, η and maximized with respect to the Lagrange multipliers α1 and α2. By carrying out this optimization, w can be written in terms of Lagrange multipliers α1 and α2. Finally, substituting the value of w and simplifying with the help of Karush-Kuhn-Tucker (KKT) the following dual problem is obtained:(3)max-12α1THHTH-1HTα1+ψ1THHTH-1HTα1-ψ1Tα1s.t.0⩽α1⩽C1e,where H=K(A,AT)e,ψ1=Y-eε1. Then, problem (3) leads to(4)-HTψ1+HTHu1+HTα1=0.That is,(5)u1=HTH-1HTψ1-α1.The solution of (5) requires the inverse of HTH. Sometimes, matrix may be ill-conditioned. This situation may be avoided by adding a regularization term, σ>0, to HTH. Here, “I” is an identity matrix of suitable dimension. Therefore, (5) is reformulated as(6)u1=w1Tb1T=HTH+σI-1HTψ1-α1.Similarly, problem (2) is considered and its dual is obtained as(7)max-12α2THHTH-1HTα2+ψ2THHTH-1HTα2-ψ2Tα2s.t.0⩽α2⩽C2e,where ψ2=Y+eε2. Hence, problem (7) leads to (8)u2=w2Tb2T=HTH+σI-1HTψ2+α2.Note that, in the duals (3) and (7), the inversion of matrix HTH of size n+1×n+1 should be computed [15]. Once the vectors u1 and u2 are known from (6) and (8) the two up and down bound functions are obtained. Then, the estimated regressor is constructed as follows:(9)ϕx=12ϕ1x+ϕ2x=12w1+w2TKA,x+12b1+b2=wmTKA,x+bm.


## 3. Identification of Wiener Models

The Wiener cascade, a linear filter followed by a static nonlinearity as shown in Figure2, is often used to represent certain higher-order nonlinear systems.Figure 2
Block diagram of Wiener model. The investigator is assumed to have access to the input,u(t), and the output, y(t), but not the intermediate signals, x(t) and y^(t), or the innovation, e(t).In this section, the development in [19] is followed, up to the point where the SVM optimization is introduced (where TSVM is used). The output error Wiener model can be described as(10)xn=Gq-1un,yn=Nxn+en,where un,xn,yn∈R are the input and output of the linear block and the measured output signals, respectively, for n=1,…,N. The innovation en is assumed to be white. Identifying output error models is computationally heavier than autoregressive exogenous (ARX) models even for linear systems. Several techniques have been suggested to deal with linear output error (OE) models such as instrumental variables (IV), subspace methods, algorithms based on the Steiglitz-McBride, and using orthonormal filter expansions like Laguerre or Kautz functions [20]. In this paper the linear part will be represented by Laguerre filters. The truncated discrete Laguerre expansion of the transfer function Gq-1 is given by(11)Gq-1=∑k=1Mgkq-11-α21-αq-1q-1-α1-αq-1k-1,for all α<1, where the function Gq-1 is assumed to be strictly proper, analytic in q>1, and continuous in q⩾1.Remark 1.
Notice that the finite impulse response FIR model can be obtained from (11) when α=0.
Substituting (11) into the first equation of (10) results in (12)xn=∑k=1Mgkq-11-α21-αq-1q-1-α1-αq-1k-1un,which can be written as(13)xn=∑k=1MgkLkun,where Lk=q-11-α2/(1-αq-1)(q-1-α)/(1-αq-1)k-1. Applying such expansion to the linear part of Wiener model enables us to rewrite the Wiener model (10) as follows:(14)yn=N∑k=1MgkLkun+en.Remark 2.
If the nonlinearity was modeled as a polynomial of orderQ(15)Nxn=∑q=1Qcqxnq,then model (14) can be seen as bank of parallel linear filters Lk’s followed by multi-input-single-output (MISO) nonlinear block as shown in Figure 3.Figure 3
Wiener model with basis function expansion of linear block.Remark 3.
Our contribution is to replace the polynomial used in (14) by TSVM which gives(16)yn=12w1+w2TKA,xn+12b1+b2+en=wmTKA,xn+bm+en,where x[n]=[x1[n],…,xM[n]]T is the output of the bank of parallel linear filters Lk’s whose components are computed using (13).Remark 4.
Model (16) might be identified by solving (3)–(8). It is clear that, solving (3)–(8), which consists of two simple quadratic programming problems without equality constraint, is easier than solving the quadratic program developed in [19].


## 4. Example 1

To show the features of the proposed algorithm, the simulation example presented in [19] is identified. Consider a valve for fluid flow control described by(17)xn-1.4138xn-1+0.6065xn-2=0.1044un-1+0.0883un-2,y0n=xn0.10+0.9xn2,where un is a pneumatic control signal applied to the valve, xn is the position of the valve plug, and y0n is the fluid flow. The fluid flow measurement yn is given by(18)yn=y0n+en,where en consists of independent sequence of normally distributed random numbers with variance 0.0025. Two input-output sequences consisting of 1000 data points were generated as follows. A pseudorandom binary signal fluctuating between −1 and +1 with a basic clock period of seven sampling intervals was produced. Then, in each time interval, the input was multiplied by a uniformly distributed random factor between 0 and 0.4. Finally, a bias 0.5 was added on it to result in a signal having an amplitude between 0.1 and 0.9. This procedure was described and used in [19]. The training and testing data sets are shown in Figures 4 and 5, respectively. The hyperparameters (M, α, ε1, ε2, C1, C2, and σ) were selected based on cross-validation method where one parameter was varied and the others were kept fixed. The value of the varying parameter that gives the least root mean square error (RMSE) value is chosen as shown in Tables 1 and 2. The cross-validation process resulted in ε1=0.04, ε2=0.04, C1=100, C2=100, and σ=1. The Laguerre filter pole and order were chosen to be α=0.1 and M=7, respectively. Figure 6 shows TSVR based algorithm estimates together with measured output for the last 100 samples of the test data. It is clear from the figure that the algorithm produced accurate result.Table 1
TSVR: root mean square errors between true and estimated output in Example 1 for test data set using various filter poles and filter orders. The other hyperparameters were fixed asε1=ε2 = 0.08, C1 = C2 = 10, and σ=1.
α
M
RMSE
CPU time
0.1
1
0.1142
19.8875
2
0.0807
19.7093
3
0.0657
18.1007
4
0.0638
22.0200
5
0.0625
23.6666
6
0.0574
18.2303
7
0.0540
18.6477
8
0.0542
22.4743
9
0.0560
21.5188
10
0.0588
25.2510
0.2
0.0548
18.6072
0.3
0.0569
24.6136
0.4
7
0.0615
19.4013
0.5
0.0833
20.5413
0.6
0.1011
18.9904Table 2
TSVR: root mean square errors between true and estimated output in Example 1 for test data set using various values of hyperparametersC1=C2,ε1=ε2, and σ. The filter pole and filter order were fixed as α=0.1, M=7, respectively.
ε
C
σ
RMSE
CPU time
0.08
10
1
0.0540
17.9305
100
0.0540
18.2326
200
0.0543
18.6875
1000
0.0539
17.9402
2000
0.0564
18.5952
0.08
100
0.1
0.1019
18.3051
0.2
0.0676
18.8480
0.04
100
1
0.0540
18.3080Figure 4
Training data of Example 1.Figure 5
Testing data of Example 1.Figure 6
Measured system output (solid line), output of identified TSVM based model (dashed line).


## 5. Example 2

Consider pH neutralization process consisting of a continuous stirred tank reactor (CSTR) where a strong base (NaOH) reacts with the feed stream, a strong acid (HCl). The process is shown in Figure7. The process input is the flow rate of the strong base and the process output is the pH of the effluent solution. The acid flow rate as well as the volume of the tank is assumed to be constant. The identification algorithm described in Section 3 is used to estimate a Wiener model of the simulation system presented in [21]. The nominal operating conditions of the system are given in Table 3. The model is highly nonlinear due to the implicit output equation, known as the titration curve. The system was excited with band limited white noise, with zero mean and 0.01 variance, around the nominal value of the base flow rate. The output of the system was corrupted with additive Gaussian white noise with zero mean and standard deviation σ= 0.001.Table 3
Nominal value.
Nominal value
Process variable
19.4 mL/s
Base flow rate (q1)
37 mL/s
Acid flow rate (q2)
0.000286 M
Base concentration
0.00015 M
Acid concentrationFigure 7
pH neutralization process.The training and testing data sets are shown in Figures8 and 9, respectively. The TSVM and SVM hyperparameters and the Laguerre filter pole and order were chosen based on cross-validation method. For TSVM, the cross-validation resulted in ε1=0.08, ε2=0.08, C1=100, C2=100, and σ=1 (see Table 5). In Table 4, the Laguerre filter pole and order were chosen to be α=0.1 and M=9, respectively. The SVM hyperparameters were chosen to be ε=0.08, C=1000, and σ=1 as shown in Table 7 and the Laguerre filter pole and order were selected as α=0.5 and M=16 (see Table 6).Table 4
TSVR: root mean square errors between true and estimated output in Example 2 for test data set using various filter poles and filter orders. The other hyperparameters were fixed asε1=ε2 = 0.08, C1 = C2 = 100, and σ=1.
α
σ
C
M
RMSE
CPU time
0.1
1
100
1
0.4466
16.3546
2
0.4125
15.5290
3
0.3831
15.8178
4
0.3465
19.4221
5
0.2821
17.2951
6
0.2439
18.0242
7
0.1953
18.2451
8
0.1462
17.7310
9
0.1248
16.7690
10
0.1272
19.6912
11
0.1311
20.2375
12
0.1313
19.3700
0.2
9
0.1275
19.6035
0.3
0.1280
20.6774
0.4
0.1287
22.1147
0.5
0.1283
22.0277
0.6
0.1303
16.1832
0.7
0.1303
21.4789
0.8
0.1334
16.3681
0.9
0.1473
18.8526Table 5
TSVR: root mean square errors between true and estimated output in Example 2 for test data set using various values of hyperparametersC1=C2, ε1=ε2, and σ. The filter pole and filter order were fixed as α=0.1, M=9, respectively.
ε
C
σ
RMSE
CPU time
0.08
1
1
0.1250
22.0614
10
0.1248
20.5972
100
0.1248
16.0475
150
0.1248
19.8612
200
0.1249
19.1683
500
0.1250
21.4734
1000
0.1248
16.5751
1500
0.1264
19.3178
2000
0.1255
20.9168
0.08
100
0.1
0.1372
21.8351
0.2
0.1353
18.4019
0.3
0.1321
20.3357
1
0.1248
16.0475
2
0.1298
23.4172
3
0.1368
21.0749
4
0.1453
19.6277Table 6
SVR: root mean square errors between true and estimated output in Example 2 for test data set using various filter poles and filter orders. The other hyperparameters were fixed asε = 0.08, C = 100, and σ=1.
α
σ
C
M
RMSE
CPU time
0.1
1
100
1
0.4533
33.3899
2
0.4291
27.5467
3
0.3944
29.6294
4
0.3655
27.5746
5
0.3074
29.2262
6
0.2651
32.6605
7
0.2318
29.8089
8
0.1926
31.2993
9
0.1498
27.3892
10
0.1483
26.7030
11
0.1472
28.4650
12
0.1466
26.9067
13
0.1455
22.8941
14
0.1444
22.9346
15
0.1439
22.8492
16
0.1433
22.1835
17
0.1435
23.0951
0.2
16
0.1429
28.7001
0.3
0.1403
28.4638
0.4
0.1371
27.4459
0.5
0.1345
28.8956
0.6
0.1364
24.4976Table 7
SVR: root mean square errors between true and estimated output in Example 2 for test data set using various values of hyperparametersC, ε, and σ. The filter pole and filter order were fixed as α=0.5, M=16, respectively.
ε
C
σ
RMSE
CPU time
0.08
1
1
0.4803
19.7045
10
0.1894
20.8680
100
0.1433
28.0950
150
0.1393
32.9487
200
0.1370
37.3656
500
0.1323
51.4939
1000
0.1311
72.5800
2000
0.1313
85.3550
0.08
1000
0.1
0.1369
73.95
0.2
0.1371
73.48
0.3
0.1371
73.19
0.4
0.1369
73.23
0.5
0.1364
73.92
1
0.1311
72.58
2
0.1643
72.46Figure 8
Training data of Example 2.Figure 9
Testing data of Example 2.By comparing the least RMSE (0.1248) and CPU time (21.36 sec.) values of the TSVR algorithm (Table5) with the RMSE (0.1311) and CPU time (72.58 sec.) values of the SVR algorithm (Table 7), it is clear that TSVM algorithm outperforms the SVM algorithm in terms of accuracy and speed of computation.


## 6. Conclusion

In this paper a new algorithm for identification of Wiener systems using TSVR has been derived and used to identify simulation examples. The algorithm was able to model the nonlinearity, without requiring any a priori assumptions regarding its structure. The linear model was represented by Laguerre filter. The main advantage of the proposed algorithm over the method presented in [19] is that the proposed algorithm results in two smaller quadratic programs which is easier to solve than the quadratic program developed in [19]. It was shown in Example 2 that TSVM algorithm is computationally lighter than SVM algorithm.



---

*Source: 125868-2015-09-16.xml*