Virtual environments already became a part of our daily life including applied computer games, learning environments (Devlin et al., 2015), business and project management environment (Mattioli et al., 2015), social networks and their games. These applications and its multimedia elements are causing negative or positive emotions and are considered as a virtual stimuli (Wrzesien et al., 2015). These stimuli may be used as a clue to regulate human psychological, emotional and social state (Devlin et al., 2015) or even to treat various mental diseases (Calvo et al., 2015). For this purpose, a control mechanism for human state regulation or stabilization is needed.

EEG-based emotion signals (excitement, frustration, engagement/boredom) are characterized as reliable and quick response signals (Hondrou and Caridakis, 2012; Mattioli et al., 2015; Sourina and Liu, 2011). However, foremost we need to compose and investigate mathematical models describing dependencies between emotion signals as a reaction to stimuli.

Predictive input-output structure models were proposed and investigated for exploring dependencies between virtual 3D face features and human reaction to them when dynamic virtual 3D face is observed without a virtual reality headset (Kaminskas et al., 2014; Kaminskas and Vidugirienė, 2016; Vaškevičius et al., 2014). Predictive models are necessary in the design of predictor-based control systems (Astrom and Wittenmark, 1997; Clarke, 1994; Kaminskas, 2007). Predictor-based control was applied to the control of human emotion signals, when a dynamic 3D face is observed without a virtual reality headset (Kaminskas et al., 2015; Kaminskas and Ščiglinskas, 2016).

This paper introduces a predictor-based control with a generalized minimum variance control quality principles which are applied to the control of human response signal, when a dynamic virtual 3D face as stimuli is observed using a virtual reality headset.

A virtual 3D face with changing distance-between-eyes was used for input and EEG-based pre-processed excitement signal of a volunteer was measured as output (Fig. 1). The output signal was recorded with Emotiv Epoc device. This device records EEG inputs from 14 channels (in accordance with the international 10–20 locations): AF3, F7, F3, FC5, T7, P7, O1, O2, P8, T8, FC6, F4, F8, AF4 (Emotiv Epoc specifications, Mattioli et al. (2015)). Values of the output signal (excitement) vary from 0 to 1. If excitement is low, the value is close to 0 and if it is high, the value is close to 1.

A dynamic stimulus was formed from a changing woman face (Kaminskas et al., 2015). A 3D face created with Autodesk MAYA was used as a “neutral” one (Fig. 2, middle). Other 3D faces were formed by changing distance-between-eyes in an extreme manner (Fig. 2, top, bottom). The transitions between normal and extreme stages were programmed. “Neutral” face has 0 value, largest distance-between-eyes corresponds to value 3 and smallest distance-between-eyes corresponds to value – 3 (Fig. 2). At first “neutral” face was shown for 5 s, then the distance-between-eyes was increased continuously and in 10 s the largest distance between eyes was reached, then 5 s of steady face was shown and after that the face came back to “normal” state in 10 s. Then “normal” face was shown for 5 s, followed by 10 s of continuous change to the face with the smallest distance between-eyes, again 5 s of steady face was shown and in the next 10 s the face came back to “normal”. The experiment was continued in the same way further using 3 s time intervals for steady face and 5 s for continuous change. Eight volunteers (four males and four females) were tested. Each volunteer was watching a changing virtual 3D face with virtual reality headset and each experiment was approximately about 100 s long. EEG-based excitement and changing distance-between-eyes signals were measured with sampling period T 0 = 0.5 s and recorded in real time.

To estimate the possible relation between human response (excitement) and virtual 3D face feature (distance-between-eyes) a cross-correlation analysis was performed. The estimates of the maximum cross-correlation function values (1) max τ | r y x [ τ ] | = max τ | R y x [ τ ] R y y [ 0 ] R x x [ 0 ] | are shown in Table 1. In Eq. (19) R y x [ τ ] is a cross-covariation function between distance-between-eyes (x) and excitement (y) signals, and R y y [ τ ], R x x [ τ ] are auto-covariation functions (Vaškevičius et al., 2014). Examples of cross-correlation functions are demonstrated in Fig. 3.

The shift of the maximum values of cross-correlation functions in relation to R y x [ 0 ] allows stating that there exists dynamic relationship between virtual 3D face feature (distance-between-eyes) and human response (excitement) to them. High cross-correlation values justify a possibility of linear dynamic relations.

Dependency between human excitement signal as a response to a virtual 3D face feature (distance-between-eyes) changes is described by linear input-output structure model (Kaminskas et al., 2014) (2) A ( z − 1 ) y t = θ 0 + B ( z − 1 ) x t + ε t , where (3) A ( z − 1 ) = 1 + ∑ i = 1 n a i z − i , B ( z − 1 ) = ∑ j = 0 m b j z − j , y t is an output (excitement), x t is an input (distance-between-eyes) signals respectively observed as y t = y ( t T 0 ) , x t = x ( t T 0 ) with sampling period T 0 , ε t corresponds to noise signal, z − 1 is the backward-shift operator ( z − 1 x t = x t − 1 ) and θ 0 is a constant value.

Eq. (2) can be expressed in the following expanded form (4) y t = θ 0 + ∑ j = 0 m b j x t − j − ∑ i = 1 n a i y t − i + ε t . Parameters (coefficients b j and a i , degrees m and n of the polynomials (3) and constant θ 0 ) of the model (2) or (4) are unknown. Parameter identification is performed in accordance with the observations obtained during the experiments with the volunteers.

The current estimates of the parameters can be obtained in the identification process from the condition (Kaminskas, 2007) (5) c ˆ t : Q ˜ t ( c ) = 1 t − n ∑ k = n + 1 t ε k | k − 1 2 ( c ) → min c ϵ Ω c , where (6) c T = [ θ 0 , b 0 , b 1 , … , b m , a 1 , a 2 , … , a n ] is a vector of the coefficients of the polynomials (3) and θ 0 , (7) ε t + 1 | t ( c ) = y t + 1 − y t + 1 | t is one-step-ahead output prediction error, (8) y t + 1 | t = θ 0 + z [ 1 − A ( z − 1 ) ] y t + B ( z − 1 ) x t + 1 is one-step-ahead output prediction model, (9) Ω c = { a i : | z i A | < 1 , i = 1 , 2 , … , n } is stability domain (unity disk) for the model (2), z i A is the roots of the polynomial (10) z i A : A ( z ) = 0 , i = 1 , … , n , A ( z ) = z n A ( z − 1 ) , z is a forward-shift operator ( z y t = y t + 1 ), T is a vector transpose sign, sign | . | denotes absolute value.

Predictive model (8) can be expressed in the form of linear regression (11) y t + 1 | t = β t T c , where (12) β t T = [ 1 , x t + 1 , x t , … , x t − m + 1 , − y t , − y t − 1 , … , − y t − n ] . Considering Eq. (7) and Eq. (11), identification criterion (13) Q t ( c ) = 1 t − n ∑ k = n + 1 t ( y k − β k − 1 T c ) 2 is a quadratic form of the vector variable c.

Accordingly, solution of the minimization problem (5) is separated into two stages. In the first stage, which is application of the least squares method, parameter estimates are calculated without evaluation of restrictions (14) c t = [ ∑ k = n + 1 t β k − 1 β k − 1 T ] − 1 [ ∑ k = n + 1 t y k β k − 1 ] . In the second stage, these estimates are projected into stability domain (9) (15) c ˆ t = Γ c t , where (16) Γ = 1 0 0 0 I b 0 0 0 γ I a , 0 < γ ⩽ 1 is a diagonal block-matrix of projection of the dimension ( m + n + 2 ) × ( m + n + 2 ), I b  and I a are correspondingly unity matrix dimension ( m + 1 ) × ( m + 1 ) and ( n × n ).

Factor γ in matrix (16) is calculated by equation (17) γ = min { 1 , γ max − γ 0 } , where γ max ‖ c t ‖ is the distance from the point 0 (origin) to the boundary of stability domain Ω c in the direction of c t , ‖ . ‖ is the Euclidean norm sign, γ 0 is a small and positive constant. When n ⩽ 2 (stability domain for the model (2) is defined by linear inequations) factor γ calculation was given (Kaminskas et al., 2014) and when n ⩽ 3 (stability domain is defined by linear and quadratic inequations) factor γ calculation was given (Kaminskas and Vidugirienė, 2016).

Estimates of the model orders ( m ˆ , n ˆ ) are defined from the following conditions (Kaminskas and Vidugirienė, 2016): (18) m ˆ = min { m ˜ } , n ˆ = min { n ˜ } , where m ˜ and n ˜ are polynomial (3) degrees when the following inequalities are correct (19) | σ ε [ m , n + 1 ] − σ ε [ m , n | σ ε [ m , n ] | ⩽ δ ε , n = 1 , 2 , … , | σ ε [ m + 1 , n ] − σ ε [ m , n | σ ε [ m , n ] | ⩽ δ ε , m = 0 , 1 , … , n , σ ϵ [ m , n ] is one-step-ahead prediction error standard deviation for model order ( m , n ), δ ϵ [ 0.01 , 0.1 ] is chosen constant value (which corresponds to a relative variation of predictions error standard deviation from 1% to 10%).

Validation of the predictive models (8) was performed for each of eight volunteers (four males and four females). For the identification of unknown parameters the first 60 to 100 observations were used for each volunteer. For evaluation of the model order and prediction accuracy all 185 observations were used. Each model is selected from twelve possible models (when n = 1 , 2 , 3, m = 0 , 1 , 2 , 3). The analysis of the experiment results showed relations between distance-between-eyes and excitement and it can be described using first order ( m ˆ = 0, n ˆ = 1) model (20) y ˆ t + 1 | t = θ ˆ 0 + b ˆ 0 x t + 1 − a ˆ 1 y t .

Prediction accuracies with predictive model (20) were correspondingly evaluated using the prediction error standard deviation, relative prediction error standard deviation and average absolute relative prediction error (Vaškevičius et al., 2014): (21) σ ε = 1 N − n ∑ t = n N − 1 ( y t + 1 − y ˆ t + 1 | t ) 2 , (22) σ ˜ ε = 1 N − n ∑ t = n N − 1 ( y t + 1 − y ˆ t + 1 | t y t + 1 ) 2 × 100 % , (23) | ε ¯ | = 1 N − n ∑ t = n N − 1 | y t + 1 − y ˆ t + 1 | t y t + 1 | × 100 % , where N = 185. Parameter estimates and a predictor accuracy measures are provided in Table 2. Figures 4 and 5 show examples of one-step-ahead prediction results when we are using model (20) for four volunteers.

The analysis of the identification results showed what relations between distance-between-eyes and excitement is described by first order ( m ˆ = 0, n ˆ = 1) model (20). The validation results show that excitement can be predicted on average with less than 8% average absolute relative prediction error. Accordingly, input-output structure model (2), (3) in the predictive form (8) can be applied to the design of prediction-based control system of human excitement signal.

A predictor-based control law is synthesized by minimizing control quality criterion Q t ( x t + 1 ) in an admissible domain Ω x (Kaminskas, 2007) (24) x t + 1 ∗ : Q t ( x t + 1 ) → min x t + 1 ϵ Ω x , (25) Q t ( x t + 1 ) = E { ( y t + 1 − y t + 1 ∗ ) 2 + q x t + 1 2 } , (26) Ω x = { x t + 1 : x min ⩽ x t + 1 ⩽ x max , | x t + 1 − x t ∗ | ⩽ δ t } , where E is an expectation operator, y t + 1 ∗ is a reference signal (reference trajectory for excitement signal), x min and x max are input signal boundaries (smallest and largest distance-between-eyes), δ t > 00$]]> are the restriction values for the change rate of the input signal, and sign | . | denotes absolute value, q ⩾ 0 are weight coefficients.

Then solving the minimization problem (24)–(26) for one-step ahead prediction model (8), the control law is described by equations: (27) x t + 1 ∗ = min { x max , x t ∗ + δ t , x ˜ t + 1 } , if x ˜ t + 1 ⩾ x t ∗ , max { x min , x t ∗ − δ t , x ˜ t + 1 } , if x ˜ t + 1 < x t ∗ , (28) B ˜ ( z − 1 ) x ˜ t + 1 = − L ( z − 1 ) y t + y t + 1 ∗ − θ 0 , (29) L ( z − 1 ) = z [ 1 − A ( z − 1 ) ] , (30) B ˜ ( z − 1 ) = λ + B ( z − 1 ) , λ = q / b 0 . If the roots of polynomial (31) B ˜ ( z ) = z m B ˜ ( z − 1 ) are in the unity disk (32) | z j B | < 1 , z j B : B ˜ ( z ) = 0 , j = 1 , … , m , then from (28)–(30) the following equation is correct (33) x ˜ t + 1 = 1 b 0 + λ { ∑ i = 1 n a i y t + 1 − i − ∑ j = 1 m b j x ˜ t + 1 − j + y t + 1 ∗ − θ 0 } . If a part or all of polynomial (31) roots do not belong to the unity disk, weight factor | λ | is increased until all roots rely in the unity disk. The scheme of a generalized minimum variance controller (27)–(30) is illustrated in Fig. 6.

When inserting the control signal, which is described by equations (28) and (30), to the model (2) we get a closed-loop system equation (34) [ B ( z − 1 ) + λ A ( z − 1 ) ] y t = B ( z − 1 ) ( y t ∗ − θ 0 ) + ε t . It is clear from equation (34), what stability of the closed-loop system is dependent of characteristic polynomial (35) D ( z ) = z d D ( z − 1 ) , D ( z − 1 ) = B ( z − 1 ) + λ A ( z − 1 ) , d = max { m , n } , roots, all the roots must be inside the unity disk (36) | z i D | ⩽ 1 , z i D : D ( z ) = 0 , i = 1 , 2 , … , d . The analysis of characteristic polynomial equation (35) allows to state what having stable model in the process of the identification (5)–(10), stability of a closed-loop system is obtained with any arrangement of roots of the polynomial B ( z − 1 ), when the weight factor  | λ | is increased.

From equation (34) we get what permanent component of output signal in stationary regime ( y t ∗ = y ∗ ) is (37) y = K p ( y ∗ − θ 0 ) , where (38) K p = B ( 1 ) B ( 1 ) + λ A ( 1 ) is a gain of the transfer function of the reference signal y t ∗ in a closed – loop (39) W p ( z − 1 ) = B ( z − 1 ) B ( z − 1 ) + λ A ( z − 1 ) .

Considering the expression (38), weight factor λ is calculated by equation (40) λ = K 0 ( 1 − K p ) K p , where (41) K 0 = B ( 1 ) A ( 1 ) is a gain of the transfer function of the input-output model (2) (42) W ( z − 1 ) = B ( z − 1 ) A ( z − 1 ) .

From equation (37) follows that sistematic control error (43) e p = y ∗ − y = ( 1 − K p ) y ∗ + K p θ 0 grows if K p is significantly lower than unit (weight factor | λ | or weight coefficient q in control criterion (25) are high). Accordingly, the gain K p is selected from an interval (44) K p ∈ [ 0.8 , 1 ] , if ( b 0 > 0 ) ∧ ( K 0 > 0 ) or ( b 0 < 0 ) ∧ ( K 0 < 0 ) 0)\wedge ({K_{0}}>0)\mathrm{or}({b_{0}}<0)\wedge ({K_{0}}<0)\]]]> or (45) K p ∈ [ 1 , 1.2 ] , if ( b 0 > 0 ) ∧ ( K 0 < 0 ) or ( b 0 < 0 ) ∧ ( K 0 > 0 ) . 0)\wedge ({K_{0}}<0)\mathrm{or}({b_{0}}<0)\wedge ({K_{0}}>0).\]]]>

When K p = 1 ( λ = 0, q = 0), we get a minimum variance control, in other cases we get a generalized minimum variance control.

Modelling experiments consisted of two phases. In the first phase a human excitement signal as a response to dynamical 3D face stimuli (testing input) were observed. According with these observations parameter estimates of the predictive model (20) were calculated using identification. In the second phase, dynamical virtual 3D face features were formed according with the control law (27) and (33) (control output). The control tasks were to maintain high excitement levels (reference signals). In this case a control efficiency can be evaluated by a relative measure (46) △ y = y ¯ c − y ¯ T y ¯ T × 100 % , where y ¯ T is an average of output y t T (excitement) as a response to testing input, and y ¯ c is an average of output y t c (excitement) as a response to control input. These measures are given in Table 3. Examples of excitement control results are shown in Fig. 7 and Fig. 8 (weight factor λ = − 0.0224 and weight coefficient q = 0.00026, when K p = 0.9 or λ = − 0.2346, q = 0.00143, when K p = 0.8).

Modelling results show that using predictor-based control with constraints a sufficiently good quality of human excitement signal control can be reached. Excitement signal level can be raised up on average to about 95% (when K p = 1, minimum variance control) and about 85%–70% (when K p = 0.9 and K p = 0.8, generalized minimum variance control) in comparison with testing input.

Control quality is influenced by a control signal variation speed which is limited by the parameter δ t of the admissible domain. This parameter allows decreasing control signal variation which is usually high in minimum variance control systems without constraints. Control signal variation decreases when a generalized minimum variance control is applied. In this case, the quality of control depends on a gain coefficient in closed-loop K p (38), whose value defines weight factor λ in (30) or weight coefficient q in control criterion (25).

Experiment planning and cross-correlation analysis results demonstrated that there is a relatively high correlation between 3D face features observed using virtual reality (distance-between-eyes) and human response (excitement) to the stimuli. The shift of the maximum values of the cross-correlations functions in relation to origin allows stating that there exists linear dynamic relationship between distance-between-eyes and excitement signals. Parameter identification method for building stable input-output structure model is proposed. Identification and validation results of one-step-ahead prediction model (8) show that excitement can be predicted on average with less than 8% average absolute relative prediction error.

Accordingly, input-output structure model (2) (3) in the predictive form (8) can be applied to the design of predictor-based control system for controlling human excitement signal as a response to a dynamic virtual 3D face. Control law is synthesized by minimizing generalized minimum variance control criterion in an admissible domain for input. Calculation method of weight factor λ in control law (27)–(30) or weight coefficient q in control criterion (25) is proposed. This method is based on admissable value of the systematic control error.

Sufficiently good control quality of excitement signal, maintained signal level is at average to about 90% (when K p = 1, minimum variance control) and about 70% (when K p = 0.8, generalized minimum variance control with high weight coefficient) higher compared to testing input, is demonstrated by modelling results. Experiment results demonstrated possibility to decrease variations of the control signal using a limited signal variation speed when decreasing constant δ t in expression (27) or using a generalized minimum variance control when increasing weight factor | λ |, which is calculated according to equation (40). However, in these cases, particularly applying minimum variance control, control quality decreases.