Microsoft Word - inner cover.docx The HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY (HJIC) formerly (until 2012) the Hungarian Journal of Industrial Chemistry The HJIC is an international periodicals that focuses on results of fundamental and applied research in the field of • Biotechnology • Chemical Engineering Science • Chemical Processes • Energetics • Environmental Chemistry • Environmental Engineering & Technology • Industrial Management • Material Science • Mechanical Engineering • Mechatronics • Process & System Engineering • Recycling in the form of original papers, reviews, short communications, and conference proceedings written in English. EDITORIAL BOARD Editor-in-chief: RÓBERT K. SZILÁGYI Department of Chemistry and Biochemistry, Montana State University, Bozeman, MT, U.S.A. Honorary Senior Editor: GÉZA HORVÁTH Department of Chemical Engineering Science, University of Pannonia, Veszprém, Hungary Associate Editors: JÁNOS ABONYI Department of Process Engineering, University of Pannonia, Veszprém, Hungary DEZSŐ BODA Department of Physical Chemistry, University of Pannonia, Veszprém, Hungary NORBERT MISKOLCZI MOL Department of Hydrocarbon and Coal Processing, University of Pannonia, Veszprém, Hungary DÓRA RIPPEL PETHŐ Department of Chemical Engineering Science, University of Pannonia, Veszprém, Hungary Editors: LÁSZLÓ BARTHA MOL Department of Hydrocarbon and Coal Processing, University of Pannonia, Veszprém, Hungary KATALIN BÉLAFI-BAKÓ Research Institute of Bioengineering, Membrane Technology and Energetics, University of Pannonia, Veszprém, Hungary PETER CZERMAK Institute of Bioprocess Engineering and Pharmaceutical Technology, Mittelhessen University of Applied Sciences, Giessen, Germany DÉNES FODOR Institute of Mechanical Engineering, University of Pannonia, Veszprém, Hungary MARIA GAVRILESCU Department of Environmental Engineering and Management, Gheorghe Asachi Technical University of Iasi, Romania LÁSZLÓ GUBICZA Research Institute of Bioengineering, Membrane Technology and Energetics, University of Pannonia, Veszprém, Hungary JENŐ HANCSÓK MOL Department of Hydrocarbon and Coal Processing, University of Pannonia, Veszprém, Hungary JIŘÍ KLEMEŠ Centre for Process Integration and Intensification, University of Pannonia, Veszprém, Hungary ZOLTÁN KOVÁCS Department of Management, University of Pannonia, Veszprém, Hungary JÁNOS KRISTÓF Department of Analytical Chemistry, University of Pannonia, Veszprém, Hungary ISTVÁN SZALAI Institute of Physics and Mechatronics, University of Pannonia, Veszprém, Hungary FERENC SZEIFERT Department of Process Engineering, University of Pannonia, Veszprém, Hungary JÁNOS SZÉPVÖLGYI Research Centre for Natural Sciences, University of Pannonia, Veszprém, Hungary IMRE TÍMÁR Institute of Mechanical Engineering, University of Pannonia, Veszprém, Hungary GYULA VATAI Department of Food Engineering, Corvinus University of Budapest, Hungary GÁBOR VERESS Federation of Technical and Scientific Societies – MTESZ Budapest, Hungary IBOLYA ZSOLDOS Department of Material Science and Technology, Széchenyi István University, Győr, Hungary EDITORIAL OFFICE: UNIVERSITY OF PANNONIA, P.O. BOX 158, VESZPRÉM, 8201 (HUNGARY) Tel.: +36 (88) 624-746, e-mail: hjic@almos.uni-pannon.hu; web: hjic.mk.uni-pannon.hu Felelős szerkesztő: Szilágyi Róbert Károly Nyelvi lektor: Matthew Chen Kiadja: Pannon Egyetem, 8200 Veszprém, Egyetem u. 10. Levélcím: H-8201 Veszprém, Postafiók 158, Tel.: (88) 624-000 Felelős kiadó: a Pannon Egyetem, Mérnöki Kar dékánja Hung. J. Ind. & Chem. 2014, Vol. 42(2), pp. 57-113 Table of Contents A Comparison of Advanced Non-Linear State Estimation Techniques for Ground Vehicles GYÖRGY MAX, BÉLA LANTOS ...................................................................................................... 57-64 Telemedical Heart Rate Measurement for Lifestyle Counselling MARIO SALAI, GERGELY TUBOLY, ISTVÁN VASSÁNYI, ISTVÁN KÓSA .......................................... 65-70 The Application of Unknown Input Estimators to Damp Load Oscillations of Overhead Cranes BÁLINT PATARTICS, BÁLINT KISS ................................................................................................. 71-77 Challenges of Enterprise Policy Compliance with Smartphone Enablement or an Alternative Solution Based on Behaviour-Based User Identification SÁNDOR DOBOS, ATTILA KOVÁCS ................................................................................................ 79-84 Measurement-Based Modelling and Simulation of a Hydrogen-Generating Dry Cell for Complex Domestic Renewable Energy Systems ATTILA GÖLLEI, PÉTER GÖRBE, ATTILA MAGYAR, LÁSZLÓ NEUKIRCHNER ................................. 85-89 Quasi-Polynomial Representation-Based Control of Mechanical Systems LÁSZLÓ NEUKIRCHNER, ATTILA MAGYAR .................................................................................... 91-95 A Network Model for Simulating the Dynamic Behaviour of an Energy Distribution Network LEVENTE BÁLINT, TIBOR DULAI, ÁGNES STARK-WERNER, GERGELY APOSTOL, ÁDÁM TOMPOS ........................................................................................................................... 97-101 On the Parametric Uncertainty of Weakly Reversible Realizations of Kinetic Systems GYÖRGY LIPTÁK, GÁBOR SZEDERKÉNYI, KATALIN M. HANGOS .............................................. 103-107 Stability and Parameter Sensitivity Analyses of an Induction Motor ATTILA FODOR, ROLAND BÁLINT, ATTILA MAGYAR, GÁBOR SZEDERKÉNYI ........................... 109-113 EDITORIAL PREFACE Systems and control theory is a constantly evolving scientific area and the dominant driving force in key industries and engineering fields e.g. process engineering, automotive engineering, bioengineering, and the energy industry. The aim of the current issue (Volume 42, Number 2) is to provide an overview of research topics pursued by selected PhD students. The papers presented here were selected from contributions at the 13th International PhD Workshop on Systems and Control Conference held on August 25, 2014 (virt.uni-pannon.hu/phdws2014). The objective of the conference was to establish an international forum for young researchers. The meeting provided opportunities for the participants to present and discuss the latest results and up-to-date applications in systems and control. This issue represents the entire spectrum of systems and control engineering as follows: • process modelling and analysis • control (traditional, intelligent, adaptive, etc.) • process monitoring and supervision • system identification and signal processing • bioengineering • traffic control • reaction kinetic networks • modelling of complex systems (classical, hierarchical, Bayesian, fuzzy, networks) • image processing and pattern recognition • artificial intelligence • soft computing (neural, genetic, fuzzy algorithms, etc.) • software (parallel computing, distributed and network computing, data visualization) • decision making (decision support, data mining) • applications of systems and control theory The organizers are grateful for the contributions of the authors. The tradition of the International PhD Workshop continues. You are invited to participate at the 14th International PhD Workshop in Veszprém in 2015! ATTILA MAGYAR University of Pannonia, Veszprém, HUNGARY Guest Editor HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 42(2) pp. 57–64 (2014) A COMPARISON OF ADVANCED NON-LINEAR STATE ESTIMATION TECHNIQUES FOR GROUND VEHICLES GYÖRGY MAX! AND BÉLA LANTOS 1 Department of Control Engineering and Informatics, Budapest University of Technology and Economics, Magyar tudósok körútja 2, Budapest, 1117, HUNGARY !Email: max@iit.bme.hu This paper presents a comparative examination of state-of-the-art non-linear state estimation techniques. First, the Sigma- Point (SPKF) stochastic estimators are investigated on the basis of the Unscented Transformation. This is followed by the technique of non-linear symmetry-preserving observers for non-linear dynamics on the basis of symmetrical LIE groups. The variety of the deterministic and stochastic state estimation algorithms is studied the kinematics of ground vehicles. Two cases are described that deal with the reliability of path tracking, the convergence and speed of off-line calibration. Keywords: state estimation, invariant observers, unscented KALMAN filter, vehicle kinematics Introduction Vehicles are indispensable in modern society, and vehicle safety is consequently of importance in everyday life. A lack of information about the state of vehicle and parameters presents a major obstacle for the development of vehicle control systems. The effectiveness of vehicle stability control largely depends on the accuracy of the vehicle state parameters that are measured by appropriate sensors. Unfortunately, measurements from these sensors contain bias, as well as electrical noise, and can drift with temperature changes. However, the state estimation method is also able to dramatically reduce errors introduced by measurement and process noise contained in the signal. The focus of this work is the non-linear state estimation techniques. Firstly, the ‘Sigma-Point KALMAN filter family’ (SPKF) is discussed on the basis of the ‘Unscented Transformation’ (UT) and STERLING’s interpolation [3, 5, 6]. Secondly, the symmetry-preserving observer is presented [15, 16, 18], which is based on Lie-transformation of invariant frames. The various state estimation procedures are then examined and compared through an illustrative non- linear example of vehicle kinematics. Sigma-Point KALMAN Filters Sigma-Point KALMAN filter [2, 4, 7] is a collective name for derivativeless KALMAN filters that employs the deterministic sampling based Sigma-point approach to calculate the optimal terms of the state estimation. Common methods use the Scaled Unscented Transfor- mation [4] and STERLING’s polynomial interpolation. The former replaces the original set of sigma-points with a transformed set in order to minimize the errors of higher order terms, the latter approximates the derivation of the divided difference filter [10, 11] and central difference filter in [12]. The Scaled Unscented Transformation The Scaled Unscented Transformation (SUT) has two steps. First, a fixed number of sigma points is chosen to capture the desired moments (at least mean and covariance) of the original distribution. Then the sigma points are propagated through the non-linear function and the moments of the transformed variables are estimated. The advantage of SUT over the Taylor series based approximation is that SUT is better at capturing the higher order moments caused by the non-linear transformation, as discussed in [1, 8-10]. Consider a non-linear function y = g(x) and assume x has mean x and covariance xP . The first two moments of y are calculated by the following two steps: 1) Weighted Sigma-point selection: X0 = x Xi = x + (L +λ)Px( )i=1…L Xi = x − (L +λ)Px( )i=L+1…2L (1) with weights w0 m = λ L+λ w0 c = λ L+λ + (1−α2 +β ) wi m = wi c = λ 2(L+λ ) (2) 58 where wi m, and wi c are associated with the ith sigma- point such that 1 2 0 =∑ = L i iw . Symbols λ, α, and β denote scaling parameters and ( ) ixPL )( λ+ is the ith column of the matrix square root of the weighted xPL )( λ+ covariance matrix. 2) Propagation: )( ii XgY = y ≈ wi mYi i=0 2L ∑ Py ≈ wi c(Yi − y )(Yi − y )T i=0 2L ∑ Pxy ≈ wi c(Xi − x )(Yi − y )T i=0 2L ∑ (3) where y , yP , and xyP are the approximated mean, covariance and cross-covariance respectively. The Unscented KALMAN Filter (UKF) The algorithm of the UKF according to Refs. [3, 7, 8] consists of the following steps: 1) Initialization: x̂0 = E x0[ ] , Px0 = E (x0 − x0 )(x0 − x0 ) T" # $ % (4) x̂0 a = E x̂0 a! " # $= x̂0 T v0 T n0 T! "% # $& T (5) P0 a = Px0 0 0 0 Rv 0 0 0 Rn ! " # # # # $ % & & & & (6) 2) State estimation for ∞= …1k : a) Calculate Sigma-points: Xk−1 a = x̂k−1 a x̂k−1 a +γ Pk−1 a x̂k−1 a −γ Pk−1 a" # $ % & ' (7) b) Time-update: Xi,k|k−1 x = f (Xi,k−1 x ,Xk−1 v ,uk−1), i = 0,…, 2L x̂k − = wi mXi,k|k−1 x i=0 2L ∑ Pxk − = wi c(Xi,k|k−1 x − x̂k −)(Xi,k|k−1 x − x̂k −)T i=0 2L ∑ (8) c) Measurement-update: Yi,k|k−1 = h(Xi,k|k−1 x ,Xk−1 n ), i = 0,…, 2L ŷk − = wi mYi,k|k−1 i=0 2L ∑ Pŷk = wi c(Yi,k|k−1 − ŷk −)(Yi,k|k−1 − ŷk −)T i=0 2L ∑ Pxkyk = wi c(Xi,k|k−1 x − x̂k −)(Yi,k|k−1 − ŷk −) i=0 2L ∑ (9) Kk = Pxkyk Pŷk −1 x̂k = x̂k − +Kk (yk − ŷk −) Pxk = Pxk − −KkPŷkKk T (10) where xk a = xk T vk T nk T! "# $ %& T and Xk a = (Xk x )T (Xk v )T (Xk n )T! "# $ %& T are the L dimension vectors of the augmented states and Sigma-points respectively, vR and nR are the covariances of the process and measurement noise respectively, λγ += L , λ = α2(L+κ)+L, α = 10-3, κ = 0, β = 2, and wi is the weighting parameter defined in Eq.(2). The Central-Difference KALMAN Filter This approach is based on STERLING’s polynomial interpolation, which was used in the derivations of divided difference filter [10, 11] and central difference filter in Ref. [12]. The result of this approximation can be used for sigma-point approach as in the case of SUT. The algorithm consists of the following steps: 1) Weighted Sigma-point selection: X0 = x Xi = x + h Px( )i=1…L Xi = x − h Px( )i=L+1…2L (11) with weights w0 m = h2−L h2 , wi m = 1 2h2 , wi c1 = 1 4h2 , wi c2 = h2−1 4h2 , (12) where wi m, and wi c are associated with i = 1, …, 2L sigma-point. 2) Initialization: the same as defined in Eqs.(4-6) 3) State estimation for k = 1, …, ∞: a) Sigma-points for time-update: x̂k−1 av = x̂k−1 v" # $ % , Pk−1 av = Pxk−1 0 0 Rv " # $ $ % & ' ' (13) Xk−1 av = x̂k−1 av x̂k−1 av + h Pk−1 av x̂k−1 av − h Pk−1 av " # $ % & ' (14) b) Time-update: Xi,k|k−1 x = f (Xi,k−1 x ,Xk−1 v ,uk−1), i = 0,…, 2L x̂k − = wi mXi,k|k−1 x i=0 2L ∑ Pxk − = wi c1 (Xi,k|k−1 x − XL+i,k|k−1 x )2 + i=1 L ∑ wi c2 (Xi,k|k−1 x + XL+i,k|k−1 x − 2X0,k|k−1 x )2 i=1 L ∑ (15) 59 c) Sigma-points for measurement-update: x̂k−1 an = x̂k − n" #$ % &' , Pk|k−1 an = Pxk − 0 0 Rn " # $ $ % & ' ' (16) Xk|k−1 an = x̂k|k−1 an x̂k|k−1 an + h Pk|k−1 an x̂k|k−1 an − h Pk|k−1 an " # $ % & ' (17) d) Measurement-update: Yi,k|k−1 = h(Xi,k|k−1 x ,Xk|k−1 n ), i = 0,…, 2L ŷk − = wi mYi,k|k−1 i=0 2L ∑ Pŷk = wi c1 (Yi,k|k−1 −YL+i,k|k−1)2 i=0 2L ∑ + wi c2 (Yi,k|k−1 +YL+i,k|k−1 − 2Y0,k|k−1)2 i=0 2L ∑ Pxkyk = w1 c1Pxk − (Y1:L,k|k−1 −YL+1:2L,k|k−1)T Kk = Pxkyk Pŷk −1 x̂k = x̂k − +Kk (yk − ŷk −) Pxk = Pxk − −KkPŷkKk T (18) Symmetry-Preserving Observers Among many problems in control theory, the symmetric nature of the system can be utilized, such as in the case of optimal control for feedback or regulations. The symmetry properties can also be used in the design of observers. Symmetry-preserving observers are also known as invariant observers that are based on the differential geometric background of abstract LIE- groups [15, 16, 18]. The constructive algorithm can be defined for designing invariant observers, which is based on the group transformation of symmetry. The aim is to produce invariant frames and invariant output errors to transform a locally asymptotically convergent observer around an equilibrium point into an invariant one, retaining its first-order approximation. The method in Refs. [16, 18] benefits from the fact that the error equations and the first-order approximation can be calculated explicitly and the relationships can be defined globally. In addition, the invariant error equations make stability issues easier to deal with [19- 21]. Assumptions Suppose G is a LIE-group with identity e on manifold Σ. The group transformation (ϕg )g∈G on Σ is a smooth mapping (g,ξ )∈G ×Σ→ϕg(ξ )∈ Σ (19) such that ϕe(ξ ) = ξ and ϕg2 (ϕg1 (ξ )) =ϕg2g1 (ξ ) for ζ, g1, and g2. The transformation group is local if ϕg(ξ ) is defined only when g lies sufficiently near to e. We consider only local transformations. Invariant observers are designed for the smooth non- linear system d dt x = f (x,u) , y = h(x,u) (20) where x, u, and y belong to the open subsets X ⊂ ℝn, ⊂U ℝm, and ⊂Y ℝp respectively and p ≤ n . Let r ≤ n be the dimension of the LIE-group. We assume that for each x the mapping x→φg(x) is full rank. Time functions u(t) and y(t) are known. Invariant Pre-Observer For the differential geometric description of the symmetric observers [14], the following definitions and theorems should be introduced: Definition 1. The group of local transformation on manifold UX× is expressed by (X,U) = (φg(x),ψg(u)) (21) where φg(x),ψg(u) is local diffeomorphism and g ∈G Definition 2. System d dt x = f (x,u) is G-invariant, if f (φg(x),ψg(u)) = Dφg(x) ⋅ f (x,u) for g, x, u. Thus, the system remains unchanged for the local transformation, !X = f (X,U) . Definition 3. The output y = h(x,u) is G-equivariant, if there exists ρg transformation on Y, such that g ∈G and h(φg(x),ψg(u)) = ρg(h(x,u)) for g, x, u. Thus, the output remains unchanged for the local transformation, Y = h(X,U) . Definition 4. The vector field w is G-invariant on manifold X, if d dt x = w(x) is invariant to the transformation, i.e. w(φg(x)) = Dφg(x) ⋅w(x) for g, x. Definition 5. An invariant frame (w1, w2, …, w3) consists of n linearly point-wise independent G- invariant vector fields on X thus (w1(x),w2(x),…,wn (x)) is a basis of the tangent space TX(x) . Lemma 1. The invariant frame can be constructed from the canonical basis of X. The vector field defined by wi(x) = (Dφγ (x) ) −1 ⋅ ∂∂xi (22) for i = 1, …, n forms an invariant frame. Definition 6. The system d dt x̂ = F(x̂,u, y) is a pre- observer, if F(x,u,h(x,u)) = f (x,u) is satisfied for x and u. Furthermore, if lim x̂(t) = x(t) as t→∞ , then the pre-observer is asymptotic. 60 Definition 7. The pre-observer is G-invariant, if ),,ˆ()ˆ())(),(),ˆ(( yuxFxDyuxF gggg ⋅= ϕρψϕ (23) is satisfied for g, x̂ , u, and y i.e. the observer remains unchanged for the group transformation on YU,,X . Thus, X = φg(x) , U =ψg(u) , Y = ρg(y) , and d dt X̂ = F(X̂,U,Y ) . An invariant observer is an asymptotic G-invariant pre-observer. Moreover, if the pre-observer is G- invariant and rank (F) = dim y , then the output map y is G-equivariant. Definition 8. Smooth map (x̂,u, y)→ E(x̂,u, y)∈ ℝp is an invariant output error if a) y→ E(x̂,u, y) is invertible for x̂ , u, and y b) E(x̂,u,h(x̂,u)) = 0 for x̂ and u c) E(φg(x̂),ψg(u),ρg(y)) = E(x̂,u, y) for x̂ , u, and y. The first two properties mean E is an output error, the third expresses the invariance. The introduction of the invariant error is necessary, because the output error ŷ− y does not usually preserve the geometry of the system [21]. The following three theorems [15, 16] support the construction of the symmetry-preserving observers. Theorem 1. The observer d dt x̂ = F(x̂,u, y) is a G- invariant pre-observer for the G-invariant non-linear Eq.(20) with a G-equivariant output if and only if F(x̂,u, y) = f (x̂,u)+ Li(I(x̂,u),E(x,u, y)) i=1 n ∑ wi(x̂) (24) where E is an invariant output error, (x̂,u)→ I(x̂,u)∈ ℝn+m&r is a full rank function, Li is a smooth function, such that Li(I(x̂,u), 0) = 0 and (w1,w2,…,wn ) is an invariant frame. Since Li(I,E) = Li(I,E)E Li(I,E)wi i=1 n ∑ = wi(Li(I,E) i=1 n ∑ E) = w1 … wn( ) L1(I,E) ! Ln (I,E) " # $ $ $ $ % & ' ' ' ' E (25) The observer can be written as follows F(x̂,u, y) = f (x̂,u)+W (x̂)L(I(x̂,u),E(x̂,u, y))E(x̂,u, y). (26) Theorem 2. The following statements are valid: • (x̂,u, y)→ E(x̂,u, y) invariant output error exists, • there is a (x̂,u)→ I(x̂,u)∈ℝn+m&r invariant function, • every other invariant output error has the form of !E(x̂,u, y) = L(I(x̂,u),E(x̂,u, y)) (27) where L is a smooth function such that L(I, 0) = 0 and E →L(I,E) is invertible. Theorem 3. Consider the equilibrium of a non-linear system f (x,u ) = 0 and y = h(x,u ) Assume that the linearized system A,B,C,D around this equilibrium is observable, where A = ∂f ∂x (x,u ) , B = ∂f ∂u (x,u ) ,C = ∂h ∂x (x,u ) , A = ∂h ∂u (x,u ) and let L be such that A+ LC is a stable matrix. From the locally asymptotically convergent observer d dt x̂ = f (x̂,u)+ L(ŷ− y) , ŷ = h(x̂,u) (28) an invariant observer can be constructed by the same linear approximation d dt x̂ = f (x̂,u)+W (x̂)L(I(x̂,u),E(x̂,u, y))E(x̂,u, y) (29) where L = −W (x)−1LV−1 and V = ∂E ∂y(x,u,y ) is an invertible matrix of pp× dimensions. There is no general algorithm for designing the Li gain functions of Theorem 1. However, if we consider the following invariant state error η(x, x̂) = φγ (x)(x̂)−φγ (x)(x) (30) instead of the xx −ˆ state error, where )(xγ is the solution of φg a (x) = c (31) with respect to g, then the invariant error will only depend on the I(x,u) trajectories. Constructive Algorithm Consider an invariant system, i.e. unchanged by transformation Eq.(21) with an equivariant output (Definition 3). The non-linear observer design for non- linear dynamic systems [10, 11] can be divided into the following steps: a) Choose a G LIE-group and a group transformation taking into account the symmetrical properties of the system. b) Solve the normalization Eq.(31). Build an invariant error E and the complete set of scalar invariants I. c) Construct the invariant frame Eq.(22). d) Determine the pre-observer Eq.(25). e) Linearize the system f (x,u ) = 0 around the equilibrium and obtain A,B,C,D matrices. Check the observability and design an invariant observer from the chosen linear observer by using Theorem 3. Choose an appropriate L . f) Choose the parameters of the linearized error equations, based on the invariant state-errors η, such that the invariant system will be asymptotically stable. 61 Then the non-linear symmetry-preserving observer is asymptotically stable, and converges locally and exponentially along all system trajectories. Case Study: A Non-Holonomic Vehicle In this section, the previously introduced state estimation procedures are employed for an example of a non-holonomic vehicle tracking control. The well- known model in the robotics of the simplified vehicle kinematics is considered. The tracking is maintained by satisfying the kinematic constraints with the appropriate choice of reference signals. Kinematic Model The simplified kinematic model of the vehicle can be written in the following form !x !y !θ ! " # # # $ % & & & = ucosθ usinθ uv ! " # # # $ % & & & , h(x, y,θ ) = (x, y) (32) where u is the velocity, v is the curvature, i.e. a function of the steering angle and h is the measured output. All signals contain additive Gaussian noise. The model assumes that the contact point of the wheel and the ground does not slip. Notice, that the linearized system is not observable. Tracking For path tracking the control signals are set to satisfy the kinematic constraints, thus u and v are calculated as follows: u(t) = !xr 2(t)+ !yr 2(t) and v(t) = !!yr (t) !xr (t)− !xr (t)!!yr (t) u3(t) (33) where xr(t) and yr(t) are time functions of the reference signal. The initial conditions for tracking control is θr (0) = tan −1 !yr (0) !xr (0) " # $ % & ' . (34) Hence, the reference signals for the vehicle velocity and steering function can be calculated. In the simulations, sinusoidal path tracking is implemented using the reference signals xr (t) = t, and yr (t) = Asin(ωt). Invariant Observer Design Group Transformation Eq.(32) is independent of the origin and of the orientation of the frame chosen, i.e. it is invariant with regards to the group of rotations and translations. We chose the Lie-group G = SE(2) and defined the group transformation as follows: φ(xg,yg,θg )(x, y,θ ) = xg yg θg ! " # # # # $ % & & & & ⋅ x y θ ! " # # # $ % & & & = = x cosθg − ysinθg + xg xsinθg + ycosθg + yg θ +θg ! " # # # # $ % & & & & (35) ψ(xg,yg,θg )(u,v) = u v ! " # $ % & ρ(xg,yg,θg )(x, y) = x cosθg − ysinθg + xg xsinθg + ycosθg + yg ! " # # $ % & & One can check that Eq.(32) is invariant to the transformations ϕ, and ψ in terms of Definition 1 and the output is equivariant to transformation ρ in terms of Definition 3. Invariant Output Errors The normalization Eq.(31) for :0=c 0 = x cosθγ − ysinθγ + xγ 0 = xsinθγ + ycosθγ + yγ 0 =θ +θγ (36) The solution of Eq.(36), the scalar invariants and output errors can be written as follows respectively xγ yγ θγ ! " # # # # $ % & & & & = x y θ ! " # # # $ % & & & −1 = −x cosθ − ysinθ xsinθ − ycosθ −θ ! " # # ## $ % & & && = γ x y θ ! " # # # $ % & & & (37) I(x, y,θ,u,v) =ψγ (x,y,θ )(u,v) = u v ! " # $ % & (38) E = ρ(xγ ,yγ ,θγ )(x̂, ŷ)− ρ(xγ ,yγ ,θγ )(x, y) = cosθ̂ sinθ̂ −sinθ̂ cosθ̂ " # $ $ % & ' ' x̂ − x ŷ− y " # $$ % & '' (39) Invariant Frame W = (Dφγ (x,y,θ )(x, y,θ )) −1 ∂ ∂(x,y,θ ) = = Dφγ (x,y,θ )(x, y,θ ) = cosθ −sinθ 0 sinθ cosθ 0 0 0 1 # $ % % % & ' ( ( ( . (40) 62 Invariant pre-observer d dt x̂ ŷ θ̂ ! " # # # $ % & & & = ucosθ̂ usinθ̂ uv ! " # # # # $ % & & & & + cosθ̂ −sinθ̂ 0 sinθ̂ cosθ̂ 0 0 0 1 ! " # # # # $ % & & & & L × cosθ̂ sinθ̂ −sinθ̂ cosθ̂ ! " # # $ % & & x̂ − x ŷ− y ! " ## $ % && (41) where L is a smooth 23× gain matrix, whose elements depend on the invariant error E and the invariants I. Error Equations and Error Dynamics η = ηx ηy ηθ ! " # # # # $ % & & & & = γ (x̂, ŷ,θ̂ ) x̂ ŷ θ̂ ! " # # # $ % & & & −γ (x̂, ŷ,θ̂ ) x y θ ! " # # # $ % & & & = = (x̂ − x)cosθ̂ + (ŷ− y)sinθ̂ −(x̂ − x)sinθ̂ + (ŷ− y)cosθ̂ (θ̂ −θ ) ! " # # # # $ % & & & & (42) !ηx !ηy !ηθ ! " # # # # $ % & & & & = u(1− cosηθ )+ (uv+ L31ηx + L32ηy )ηy usinηθ − (uv+ L31ηx + L32ηy )ηx 0 ! " # # # # $ % & & & & + +L ηx ηy ! " # # $ % & & (43) It can be seen that the invariant error equations are independent of the trajectory and only depend on the relative quantities ηx, ηy, and ηθ and the I invariants. Stability and Convergence The weighting matrix is defined as follows L = − |u | a ubEy −uv uv−ubEy − |u | c 0 −ub " # $ $ $ $ % & ' ' ' ' , (44) where a and b are care positive scalars and Ey = ηy. The invariant error dynamics can be made locally asymptotically stable. Furthermore, the resulting symmetry-preserving observer is almost globally asymptotically convergent, i.e. it converges any initial condition except for one as described in Refs. [15-16]. Simulation Results The state estimation techniques are illustrated in two tracking examples. For comparative purposes, the following algorithms have been implemented: Extended KALMAN Filter (EKF), Central Difference and Unscented KALMAN Filter (CKF, UKF) and the symmetry preserving observer (SYM) discretized by the sample time of the sigma-point estimators. The reference paths are xr (t) = t, and yr (t) = Asin(ωt) , where A = 3 m and ω = π/5 rad s-1. Zero mean additive Gaussian noise is assumed with standard deviations of 0.01 and 0.005 for the process noise and 0.02 for the measurement noise. The sampling time is Ts = 0.1 s. The initial values are x0 est = y0 est =θ0 est =1 , and x0 = y0 = 0 . Example 1 The vehicle starts from the origin of the xy-plane with large initial estimation error and tracks a sinusoidal path. Fig.1 shows the estimated output positions of the vehicle using sigma-point KALMAN filters and symmetry-preserving observers. The estimation errors of different algorithms can be seen in Figs.2 and 3. The errors are smaller and show smoother signals in the case of the SPKF estimators, which was expected, since the SYM observer is deterministic. Large errors appear at the peak values of the reference signals; however, both estimation techniques achieved the same satisfactory results in spite of the relatively small sampling frequency. Figure 2: Estimation errors of the sigma-point and extended KALMAN filters from Example 1 63 Root-mean-square errors can be seen in Table 1. We can conclude that among the processes of state estimation techniques the SYM observer is competitive with SPKF estimators. Example 2 In this example the vehicle performs an off-line estimation before movement starts, followed by an on- line estimation along a sinusoidal path. It is notable that in the stationary state the velocity is zero and at the beginning of movement the velocity changes to a non- zero value abruptly, see Eq.(33). Figs.4-7 illustrate the estimated output positions and the measurement errors for the sigma-point estimators and the invariant observers respectively The error signals show a quick convergence in the case of the symmetry-preserving observer. The speed of convergence depends on the parameters set in the L matrix; however, high gains can cause the observer to be divergent and make larger errors. Since there are no u and δ signals in the stationary position and the orientation is not directly measured, the orientation can only be estimated after the vehicle starts moving. The error signal oscillates around zero in a similar way to Example 1. Due to calibration, smaller errors appear. The invariant observer and sigma-point estimators have the same magnitude of error components. The root-mean-square errors of the simulation can be seen in Table 2. Figure 3: Estimation errors of the symmetry-preserving observer (SYM) from Example 1 Figure 6: Estimation errors of the sigma-point and the extended KALMAN filters from Example 2 Figure 4: Estimated output position of the vehicle using sigma-point state estimation (UKF, CKF) and extended KALMAN Filter (EKF) from Example 2 Figure 5: Estimated output position of the vehicle using symmetry-preserving observer (SYM) from Example 2 Figure 7: Estimation errors of the symmetry- preserving observer (SYM) from Example 2 Table 1: Root-mean-squared errors of the estimated signals from Example 1 RMSE x, m y, m θ, rad UKF 0.1368 0.1587 0.0338 CKF 0.1366 0.1577 0.0336 EKF 0.1371 0.1551 0.0352 SYM 0.1371 0.1628 0.0455 Table 2: Root-mean-square errors of the estimated signals from Example 2 RMSE x, m y, m θ, rad UKF 0.0998 0.1180 0.0627 CKF 0.0993 0.1175 0.0628 EKF 0.1013 0.1165 0.0648 SYM 0.1264 0.1453 0.0760 64 Conclusion In this paper we presented a comparison of two different recently developed state estimation techniques. The Sigma-Point (SPKF) stochastic estimator is based on the Unscented Transformation and the non-linear symmetry-preserving observers are based on the group transformation of symmetrical LIE groups and invariant frames. Two examples are considered for the estimation of vehicle kinematics: the first example deals with the path tracking reliability of the vehicle; the second example examines the convergence and the speed of the off-line calibration followed by path tracking. The examples illustrate that for non-linear systems in noisy environments the use of stochastic state- estimators is needed. Moreover the Sigma-point approach provides more accurate estimation results than the commonly used non-linear EKF state estimators. The symmetry preserving observers for non-linear systems benefit from their constructive design method, which provides local asymptotic convergence. This method can be utilized by non-linear applications in which the physical and mechanical symmetry properties of the system can be exploited. REFERENCES [1] BAR-SHALOM Y., LI X.-R., KIRUBARAJAN T.: Estimation with Applications to Tracking and Navigationm, Wiley Interscience, New York, NY, U.S.A., 2001 [2] SÄRKKÄ S.: Recursive Bayesian inference on stochastic differential equations, PhD Thesis, Helsinki University of Technology, Helsinki, Finland, 2006 [3] JULIER S.J., UHLMANN J.K.: Unscented filtering and non- linear estimation, Proc. IEEE, 2004, 92(3), 401–422 [4] VAN DER MERWE R., WAN E.A., JULIER S.: Sigma-Point KALMAN Filters for Nonlinear Estimation and Sensor- Fusion: Applications to Integrated Navigation, Proc. AIAA Guidance, Navigation & Control Conference (GNC), Providence, RI, U.S.A. 2004 [5] VAN DER MERWE R.: Sigma-Point KALMAN Filters for Probabilistic Inference in Dynamic State-Space Models, PhD Thesis, Oregon Health & Science University, Eugene, OR, U.S.A., 2004 [6] JULIER S.J., UHLMANN J.K.: Corrections to unscented filtering and non-linear estimation, Proc. IEEE, 2004, 92(12), 1958–1958 [7] JULIER S.J., UHLMANN J.K.: Unscented filtering and non- linear estimation, Proc. IEEE, 2004, 92(3), 401–422 [8] WAN E.A., VAN DER MERWE R.: The unscented KALMAN filter, in KALMAN Filtering and Neural Networks, Ed.: HAYKIN S., John Wiley and Sons, New York, NY, U.S.A., 2001, Chapter 7 [9] NORGAARD M., POULSEN N., RAVN O.: Advances in Derivative-Free State Estimation for Non-linear Systems Tech. Rep. IMM-REP-1998-15, Dept. Math. Modelling, Technical University of Denmark, Lyngby, Denmark, 2000 [10] NORGAARD M., POULSEN N., RAVN O.: New Developments in State Estimation for Nonlinear Systems, Automatica 2000, 36(11), 1627–1638 [11] ITO K., XIONG K.: Gaussian Filters for Nonlinear Filtering Problems, IEEE Trans. Automatic Control 2000, 45(5), 910–927 [12] VAN DER MERWE R., WAN E.: Efficient Derivative-Free KALMAN Filters for Online Learning, Proc. 9th Eur. Symp. Artificial Neural Networks (ESANN) Bruges, Belgium, 2001, 205–210 [13] VAN DER MERWE R., WAN E.: The Square-Root Unscented KALMAN Filter for State- and Parameter- Estimation, Proc. IEEE Int. Conf. Acoustics, Speech, and Signal Processing (ICASSP), Salt Lake City, UT, 2001, 6, 3461–3464 [14] SALAM ERAN: Filtering Algorithms and Avionics Systems for Unmanned Aerial Vehicles, PhD Thesis, MINES ParisTech, Paris, France, 2010 [15] BONNABEL S., MARTIN P., ROUCHON P.: Non-linear Symmetry Preserving Observers on LIE-Groups, IEEE Transactions on Automatic Control, 2009, 54(7), 1709- 1713 [16] BONNABEL S., MARTIN P., ROUCHON P.: Symmetry- preserving Observers, IEEE Transactions on Automatic Control, 2008, 53(11), 2514–2526 [17] BONNABEL S., MARTIN P., ROUCHON P.: A non-linear symmetry-preserving observer for velocity aided inertial navigation, Proc. 2006 American Control Conference, Minneapolis, MN, U.S.A., 2006 [18] AYHANMAN N., ROUCHON P.: On invariant asymptotic observers, CDC, 2002 [19] ROUCHON P., GUILLAUME D.: Observation and Control of Simplified Car, IFAC, 1998 [20] BONNABEL S.: Symmetries in observer design: review of some recent results and application to EKF-based SLAM, Robot Motion and Control 2011, Lecture Notes in Control and Information Sciences 2012, 422, 3–15 [21] COLLON C., RUDOLPH J., WOITTENNEK F.: Invariant Feedback Design for Control Systems with Lie- symmetries – a kinematic car example, Discrete Cont. Dyn. Sys. 2011, 2(11), 312–321 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 42(2) pp. 65–70 (2014) TELEMEDICAL HEART RATE MEASUREMENTS FOR LIFESTYLE COUNSELLING MARIO SALAI,! GERGELY TUBOLY, ISTVÁN VASSÁNYI, AND ISTVÁN KÓSA Faculty of Information Technology, University of Pannonia, Egyetem u. 10.,Veszprém, 8200, HUNGARY !Email: mario.salai@gmail.com In this paper we analyse a low-cost commercial chest belt to be integrated into a lifestyle counselling system as a source of heart rate data. We compared data from a Schiller ECG Holter device, which serves as a reference to a CardioSport device. Due to missing data in the CardioSport device caused by loss of contact with the body, the creation of special algorithms was necessary for synchronization and data validation. The results show that when using our synchronization algorithms the average absolute percentage error between the two signals was 2% with correlation of more than 99%. Using a data validation algorithm, we were able to get on average more than 70% of the signal with an absolute percentage error of 3% and a high average correlation of 99%. The mean RR interval values and standard deviation of RR intervals are very close to those of the reference device using both the synchronization and data validation algorithms. When using the data validation algorithm, the reference measurements produced only slightly better results with regard to false detections of atrial defibrillation than the CardioSport device. In conclusion, we found that with a simple pre- processing algorithm, CardioSport as a low-cost device can be safely integrated into a lifestyle support system as a telemedical solution. Keywords: telemedicine, lifestyle counselling, heart rate monitor Introduction Low-cost telemedical sensors are often used in modern ambient assisted living (AAL) telemonitoring and self- management systems for providing inputs to medical intelligence algorithms [1]. Such systems extend the scope of traditional health care that is based purely on data measurement. However, the proper interpretation and reliability of the results depends on the reliability of the measured data and the sensor itself. Nevertheless, there are still surprisingly few reviews reported in the literature to date on the validation of the information content of such low-cost sensors compared to the clinically accepted reference device. An example of a device that was tested for validity is the SenseWear HR Armband [2]. In this study, they used the reference device simultaneously with the tested device as a way of validating data. However, most of the compared devices are expensive high-end devices, which present an obstacle for their wide use in telemedicine. In this proof-of-concept paper, we analyse a simple commercial chest belt chosen to be integrated into the Lavinia lifestyle mirror system [3] as a source of heart rate (HR) data. In the Lavinia system, the HR signal of the patient will be used to (i) estimate the calories burnt by physical activity, (ii) calculate the heart rate variability (HRV) in order to detect periods of mental or emotional stress, and (iii) analyse arrhythmia patterns (Poincare plots) for atrial fibrillation detection. Our approach involves the comparison of the HRV and Poincare plots computed from the filtered chest belt signal, with those parameters computed from a reference Holter device. Methods Measurements Two devices were used simultaneously by a healthy volunteer over a 24 hour period. A Schiller MT- 101/MT-200 Holter device was our reference device designed for clinical use. The chest belt was a CardioSport TP3 Heart Rate Transmitter device. Since this device does not have its own memory for storing data, we used a Nexus 7 tablet with Android version 4.4.2 to connect the device via the Bluetooth 4.0 protocol and store the measured data on the tablet. Although both devices were worn by volunteers for 24 hours, only 12 hours of the overall signal were used for analysis due to frequent detachments of the device from the body during nighttime. The measurements of 12 hours were repeated on 4 additional healthy male subjects. Signal Analysis The direct comparison of measured data was not possible due to the different designs of the reference and the telemedical devices. However, we wanted to 66 compare signals directly in terms of time and also to develop a data validation algorithm for removing the noisy parts of the CardioSport device measurements reliably without using the reference data. The problem was that the chest belt was not firmly attached to the body and sudden movements of the device caused signal loss. Therefore, we needed to create a software module for synchronization and data validation before any analysis. Data validation means removing obviously bad data (artefacts) and keeping only ‘good’ data segments of sufficient length, because, as a rule of thumb, both HRV and Poincare plot computations require data chunks of at least 5 minutes. Even though the data validation algorithm removed a considerable amount of data from the original signal, we still had enough useful data for analysis from the daytime. The Synchronization Algorithm Our simple algorithm for signal synchronization uses a sliding window that passes from the beginning of the chest belt signal to the end and calculates the absolute error between the two signals. When sliding finishes, the location of the sliding window with the minimum absolute error is considered as the point where the two signals should be synchronized. This applies only if the correlation of the data in the sliding window and the same amount of data from the reference device are higher than a minimum set by the user. If these conditions are met, the algorithm copies data from the sliding window into a newly generated third signal, which represents the chest belt signal fully synchronized with the reference signal. If conditions are not met, the third signal is filled with zeros. Finally, the algorithm extracts all the highly correlated segments from the third signal ignoring zero values. Also, a file with all the merged segments is generated for general analysis. The algorithm uses the following 5 main parameters that can be set up by the user: 1. window size: amount of data copied from the signal into the sliding window (default: 200), 2. window shift step: the number of samples by which we shift the sliding window in each iteration (default: 50), 3. absolute error window: amount of data used for calculating the minimum absolute error (default: 200), 4. maximum error distance: the number of samples by which we shift the absolute error window in order to find the minimum absolute error (default: 1000), 5. minimum correlation: minimum correlation, expressed as a percentage, required for the two signals to consider data in the chest belt signal as accurate (default: 97%). Each parameter’s default value was determined empirically. After running the synchronization process, we obtained segments of highly correlated data. Fig.1 shows the distribution of the lengths of signal segments. We can see that most segments are 3 to 18 minutes long. The longest highly correlated segment with the reference data is 110 minutes long. The default parameter settings minimize the number of overly short (< 5 min) segments. Most of the bad segments (Fig.2) are shorter than one minute, and only one bad segment was 60 minutes long. Data Validation Algorithm Another type of algorithm was used in the real telemedical scenario for finding good parts of the signal without relying on reference data. This implies finding gaps and abnormal values and omitting them. First, we compared the timestamp of each data point with the timestamp of the previous one. If the difference between the timestamps was longer than 3 seconds, we marked this as a ‘gap’. The 3-second gap detection was enabled by the chest belt’s buffering system that can tolerate short detachments of the device from the body. In the second step we identified abnormal values in the signal that were treated as gaps. The abnormal values are identified by observing the mean value of 20 neighbouring data points (10 before and 10 after a given point). If the mean value differs from the value of the Figure 1: The distribution of strongly correlated segment lengths for all subjects Figure 2: The distribution of weakly correlated segment lengths for all subjects 67 current sample by more than 300 units, we consider it invalid and mark it as a gap/error in the signal. Finally, we extract the good segments from the signal with a length of more than 5 minutes. We implemented the above algorithm in a simple software tool (Fig.3). On the left-hand side we can load the two signals and set the parameter values as well as the amount of data to be analysed. The graph shows two signals after the synchronization process was completed. The user can examine signals by clicking the Previous and Next buttons. General statistics are shown in the middle part of the screen, while in the lower part, we can see the histogram, and save the histogram and results as a file. Two tabs in the top left-hand corner allow the user to switch between synchronization and data validation algorithms. Statistical Analysis Time and frequency domain analyses, correlation comparisons, mean absolute percentage errors, and the slopes of scatter plot diagrams were compared between two measurements for HRV analysis. The specificities of a self-developed atrial fibrillation detector algorithm were compared for atrial fibrillation analysis. The latter algorithm is based on the k-means clustering of Poincaré plots (consisting of RR intervals) The time and frequency domain analyses for HRV were performed using Kubios HRV analysis software, while the rest of the analysis for HRV and atrial fibrillation was performed in Microsoft Excel. Atrial fibrillation detection was done using the MATLAB environment and the results were saved as Microsoft Excel workbooks. Results and Analysis Heart Rate Variability After the synchronization process, we got strongly correlated (greater than 97%) synchronized data segments of various durations. Table 1 summarizes the duration of signals analysed. Table 2 shows results in the time domain for Schiller and CardioSport devices after using the algorithm for the synchronization of signals. Time domain analysis shows similar values for mean RR values and standard deviation (STD RR in Eq.(1)). The average Mean RR values for the Schiller and CardioSport devices are 851 and 871 respectively. The average STD RR for the Schiller device is 108 and 110 for the CardioSport device. Figure 3: Synchronization and data validation software Table 1: Signal duration after the synchronization process Subject Duration (h:m:s) #1 10:53:28 #2 8:45:40 #3 10:30:17 #4 7:46:56 68 !"#!!! = ! !!! (!!! − !!)!! !!! (1) The frequency domain analysis for the synchronization process is presented in Table 3. The absolute power was compared for very low frequencies (VLF: 0-0.04 Hz), low frequencies (LF: 0.04-0.15 Hz), high frequencies (HF: 0.15-0.4 Hz) and ratios between low frequencies and high frequencies (LF/HF). Results show no significant difference between Schiller and CardioSport device values. The average mean absolute percentage error (MAPE) between two signals is 2% with a high average correlation of close to 100%. Using the data validation algorithm, we extracted data points from the collected signals. The duration of the resulting signal is shown in Table 4. It is important to note that due to the noise on Schiller device recordings, we had to remove noisy parts from the original signal. Therefore, even though the signal was recorded continuously for 12 hours, overall duration is much less. Calculations show that in the worst scenario only 45% of the signal can be used for analysis using this data validation method, while in the best scenario this number reaches 95%. This leads to a conclusion that results are rather subject dependent. The results of data analysis in the time domain after the removal of bad parts using the validation algorithm can be seen in Table 5. The mean RR intervals for Schiller and CardioSport devices are 851 and 871 and standard deviations are 104 and 106, respectively. The CardioSport device has slightly greater values, but these are practically identical. The frequency domain analysis for the data validation process is presented in Table 6. The absolute power was compared for very low frequencies (VLF: 0- 0.04 Hz), low frequencies (LF: 0.04-0.15 Hz), high frequencies (HF: 0.15-0.4 Hz), and ratios between low frequencies and high frequencies (LF/HF). As in the synchronization process, the results show no significant difference between the Schiller and CardioSport device values. The minimum, maximum and average percentage errors on whole signals were calculated using 5 minute long sliding windows with one minute long shift steps (Table 7). Only one subject had a high maximum error value of 34%. By visual examination, it was determined that the cause of such a high error was the artefact of the Schiller device. In spite of that, the average error remained low (2%). Table 2: Time domain analysis after synchronization Mean RR (ms) a STD RR (ms) b Subject Schiller CardioSport Schiller CardioSport #1 738 755 123 125 #2 704 720 91 93 #3 908 929 90 93 #4 855 875 145 148 #5 937 959 107 109 Average 851 871 108 111 a with 2% error, b with 1-3% error Table 3: Frequency domain analysis after synchronization Schiller CardioSport Error Subject Absolute power (ms2) Absolute power (ms2) % % % % VLF LF HF LF/HF VLF LF HF LF/HF VLF LF HF LF/HF #1 7937.6 3086 1578 1.956 8444 3224 1330 2.4235 6 4 19 19 #2 5431.5 626.6 245 2.557 5723 659.3 250.9 2.6281 5 5 2 3 #3 4251.2 1927 494.4 3.898 4543 2055 538.8 3.8146 6 6 8 2 #4 12682 1790 636.5 2.813 13514 1869 621.5 3.0077 6 4 2 6 #5 6139.8 1212 476.7 2.542 6465 1274 481.4 2.6459 5 5 1 4 Table 4: Signal duration after data validation Subject Duration (h:m:s) #1 1:28:10 #2 11:20:03 #3 6:15:38 #4 9:27:07 #5 4:29:44 Table 5: Time-domain analysis after data validation Mean RR (ms) a STD RR (ms) b Subject Schiller CardioSport Schiller CardioSport #1 701 724 136 139 #2 700 717 91 93 #3 899 921 100 100 #4 846 866 139 142 #5 958 981 88 90 Average 851 871 105 106 a with 2% error, b with 0-2% error Table 6: Frequency domain analysis after data validation Schiller CardioSport Error Subject Absolute power (ms2) Absolute power (ms2) % % % % VLF LF HF LF/HF VLF LF HF LF/HF VLF LF HF LF/HF #1 10414 2297 1171 1.96 10847 2442 1004 2.43 4 6 17 19 #2 5446 631 245 2.57 5718 654 245 2.67 5 3 1 4 #3 5163 1990 523 3.80 5424 2054 540 3.80 5 3 3 0 #4 11683 1769 616 2.87 12149 1831 594 3.08 4 3 4 7 #5 4356 1235 317 3.89 4522 1303 330 3.95 4 5 4 1 69 Fig.4 represents a typical relationship between CardioSport and Schiller devices. All gradient values are close to 1. The lowest slope value is 0.98 while the highest value is 1.02. The average mean absolute percentage error (MAPE) between two signals was 3% with a strong average correlation of 99%. Atrial Fibrillation We carried out the detection of atrial fibrillation (AFib) by analysing POINCARÉ plots consisting of 30 RR intervals. We considered 30 RR intervals per iteration and in each iteration after constructing the POINCARÉ plot we calculated the dispersion around the diagonal line and used k-means based cluster analysis to determine the number of the clusters. If the dispersion was too high (greater than 0.06) and the number of clusters was 1, or the number of clusters was more than 9; we assigned “AFib” to that series of RR intervals, otherwise to “Non-AFib”. The details of the algorithm can be seen in our previous study [16]. Since our data set did not contain real AFib cases, only specificity could be calculated with regard to the efficiency of detection. The evaluation of atrial fibrillation detection results for synchronized data validation can be seen in Tables 8 and 9. Conclusion Even though the CardioSport device may suffer from signal loss due to its design, we managed to determine that it can be safely used for telemedical purposes of measuring HRV and atrial fibrillation. We found only a few usable data segments that were less than 5 minutes long. With our algorithm that detects gaps and errors in Table 7: The minimum, maximum and average percentage errors Subject Minimum error Maximum error Average error #1 0.1% 3.5% 1.5% #2 0.0% 7.7% 2.1% #3 0.0% 33.9% 3.2% #4 0.1% 6.7% 1.9% #5 0.1% 5.1% 2.2% Average 0.1% 13.4% 2.4% Figure 4: Comparison of the CardioSport and Schiller devices after data validation Table 8: Results from the synchronized data related to atrial fibrillation detection Subject Number of iterations Schiller MT-101/MT-200 system CardioSport TP3 Heart Rate Transmitter AFib cases Non-Afib cases AFib cases Non-Afib cases #1 331 26 8% 305 92% 31 9% 300 91% #2 1796 9 1% 1787 99% 3 >1% 1793 ~100% #3 1120 7 1% 1113 99% 5 1% 1115 99% #4 1427 11 1% 1416 99% 16 1% 1411 99% #5 964 46 5% 918 95% 45 5 919 95% Min - - 1% - 92% - >1% - 91% Max - - 8% - 99% - 9% - ~100% Mean - - 3% - 97% - 3% - 97% STD - - 3% - 3% - 4% - 4% Table 9: Results from the data validation process related to atrial fibrillation detection Patient Number of iterations Schiller MT-101/MT-200 system CardioSport TP3 Heart Rate Transmitter AFib cases Non-Afib cases AFib cases Non-Afib cases #1 241 3 1% 238 99% 8 3% 233 97% #2 1879 29 2% 1850 98% 2 >1% 1877 ~100% #3 808 15 2% 793 98% 3 >1% 805 ~100% #4 1296 10 1% 1286 99% 20 2% 1276 98% #5 544 6 1% 538 99% 7 1% 537 99% Min - - 1% - 99% - >1% - 97% Max - - 2% - 99% - 3% - ~100% Mean - - 1% - 99% - 1% - 99% STD - - >1% - >1% - 1% - 11% 70 the signal and removes them with an average effectiveness of more than 70%, which translates into having enough data to calculate HRV and atrial fibrillation from daytime measurements. Regarding atrial fibrillation detection, we can conclude that by using the developed data validation algorithm the reference Schiller MT-101/MT-200 measurements produced only slightly better results with regard to false detections than the CardioSport TP3 Heart Rate Transmitter. In two cases the CardioSport measurements proved to be even better than Schiller records, which implies that some relatively simple heart rate recorders are equivalent to some Holter devices after signal processing using the data validation algorithm. We have to emphasize; however, that we have not performed any measurements on actual atrial fibrillating patients yet. Therefore, the investigation of the sensitivity of our atrial fibrillation detection algorithm under the presented circumstances could be the subject of further studies. In summary, the CardioSport as a low-cost device can easily be integrated into a lifestyle support system as a telemedical solution. ACKNOWLEDGMENTS This research was supported by the European Union and co-funded by the European Social Fund “Telemedicine- focused research activities in the field of Mathematics, Informatics and Medical Sciences” TÁMOP-4.2.2.A- 11/1/KONV-2012-0073. REFERENCES [1] PATEL S., PARK H., BONATO P., CHAN L., RODGERS M.: A review of wearable sensors and systems with application in rehabilitation, J. Neuroengng. Rehab., 2012, 9(21), 1–17 [2] KRISTIANSEN J., KORSHØJ M., SKOTTE J.H., JESPERSEN T., SØGAARD K., MORTENSEN O.S., HOLTERMANN A.: Comparison of two systems for long-term heart rate variability monitoring in free- living conditions - a pilot study, Biomed Eng. Online, 2011, 10(27), 1–14 [3] KÓSA I., VASSÁNYI I., NEMES M., KÁLMÁNNÉ K.H., PINTÉR B., KOHUT L.: A fast, Android based dietary logging application to support the lifestyle change of cardio-metabolic patients, Global Telemedicine and eHealth Updates: Knowledge Resources, Eds.: M. JORDANOVA, F. LIEVENS, 2014, 7, 553–556 [4] SALAI M., TUBOLY G., VASSANYI I., KOSA I.: Reliability of telemedical Heart Rate meters, IME J., 2014, 8(5), 49–55 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 42(2) pp. 71–77 (2014) THE APPLICATION OF UNKNOWN INPUT ESTIMATORS TO DAMP LOAD OSCILLATIONS OF OVERHEAD CRANES BÁLINT PATARTICSB AND BÁLINT KISS Department of Control Engineering and Information Technology, Budapest University of Technology and Economics, Magyar Tudósok krt. 2, Budapest, 1111, HUNGARY BE-mail: patarticsbalint@gmail.com This paper focuses on the development of state estimation methods for mechanical systems with uncertain frictional pa- rameters. The goal of the study is to provide reliable angle estimation for state-feedback-based crane control solutions, designed to reduce load sway. Cranes are underactuated systems, usually unequipped with the sensors necessary to measure the swinging angle, therefore the damping of their oscillatory behaviour is a challenging task. Two estimators are proposed for the calculation of the unmeasured states. One is based on an ’unknown input Kalman filter’ (UIKF), the other applies the ’unscented Kalman filter’ (UKF) with load prediction. Simulation results are provided to demonstrate the accuracy of the algorithms. Keywords: overhead crane, state estimation, nonlinear systems, unknown inputs, Kalman filter Introduction The sway of the load carried by cranes is an unwanted phenomenon. Friction may decouple load motion from the remaining crane mechanism. This causes difficulty for inexperienced operators, especially when they try to land the cargo. The friction present in the crane mech- anism makes the reduction of the swinging particularly demanding. Therefore the goal of the design of crane con- trol systems is often to ensure minimal swinging along a specified trajectory [1–6]. There are numerous algorithms available in the liter- ature for the solution to this problem. The design of such a control system can be addressed using soft computing methods such as fuzzy logic and neural network with genetic algorithm. These solutions require however the measurement of the load position, which is often based on workspace visualization [1]. Hard commuting techniques are also widely applied for sway reduction. An H2 opti- mal solution is described in Ref.[2], which also relies on image processing for the measurement of the load coordi- nates. There are a variety of algorithms taking advantage of the flatness property (or exact linearizability) of over- head cranes [7]. The control systems in Refs.[3] and [4] consider all the states measurable, however an estimator is proposed in Ref.[3] for the computation of the uncer- tain parameters. In the case of navy cranes, local asymptotic stability can be achieved without the measurement of the swing- ing angle, if trajectory planning exploits the flatness prop- erty of the system [6]. There is also a tracking control algorithm specifically designed for overhead cranes, that is capable of effectively reducing the sway despite the friction [5]. For the calculation of the unmeasured states, this controller uses a linear observer. In this paper we propose new methods of state observation for systems that take into account friction, and are applicable to over- head cranes. These estimators can replace the linear ob- server in Ref.[5] or the need to measure the load coor- dinates in any other state-feedback-based method to en- hance the precision of the control. Our estimation meth- ods are based on the assumption that the effect of the fric- tion on the inputs can be reduced. This is true for most mechanical systems (including overhead cranes), where the actuating signals are usually forces and torques. In the following section, the modelling of Lagrangian mechanical systems is overviewed, and it is shown how the friction can be handled as an input disturbance. In the next section the state estimation techniques are presented based on the possible approaches of interpreting the fric- tion. The next section describes the results of the applica- tion of these methods to the state estimation problem of overhead cranes and presents the simulation results. Lagrangian-Based Models of Mechanical Systems The nonlinear dynamics of the controlled system can be obtained using the Euler-Lagrange equations. If the gen- eralized coordinates are q = � q1 q2 . . . qn �T and the generalized forces are ⌧ = � ⌧1 ⌧2 . . . ⌧n �T , the Euler- Lagrange equations read d dt @L @q̇i � @L @qi = ⌧i, (1) 72 where L is the so-called Lagrangian, the difference be- tween the kinetic and potential energy of the system, q̇i = dqi dt and i = 1, 2, . . . n. Using these equations, under some conditions the controlled system’s model can be written in the form H(q)q̈ + h(q, q̇) + hs(s) = ⌧, (2) where H(q) is the inertia matrix, h(q, q̇) comprises the centrifugal, Coriolis and gravitational terms, and hs(s) is the frictional term with s being the vector of the fric- tional forces. In this form hs(s) can be reorganized to the right-hand-side of the equation. Introducing the dis- turbance torques as ⌧̃ = ⌧ � hs(s), Eq.(2) yields H(q)q̈ + h(q, q̇) = ⌧̃ . (3) The state-space representation of the model is ẋ = f(x, u) y = g(x), (4) where x is the state, u is the input, y is the output of the system, and f and g are nonlinear functions. This form can be obtained from Eq.(5) assuming that xT =� qT q̇T � , and the output is not directly dependent on the torques. The state-space equations read ẋ = ✓ q̇ �H�1 (q)h(q, q̇) ◆ + ✓ 0 H�1 (q) ◆ ⌧̃ y = g(x). (5) The inputs of mechanical systems are usually forces and torques, consequently they are only dependent on ⌧̃ . This means, that the effect of the friction can be reduced to the actuated degrees of freedom of the system. If it is necessary for the design of the control system, the nonlinear dynamics can be approximated by a first order Taylor series expansion around a setpoint (x0, u0). The linear model then becomes �ẋ = A�x+B�u �y = C�x, (6) where the quantities prefixed with � mean the distance from the setpoint, e.g. �x = x� x0, and the matrices are Jacobians [8]. State Estimation Methods Since the friction can be interpreted in our setup as an input disturbance, we consider two ways for its compen- sation in the state estimation. Either all input variables can be disregarded, or the additive disturbance can be es- timated using a load predictor. These ideas form the basis of the two estimator design techniques described as fol- lows. The computations of the controllers implemented in embedded systems are based on sampled signals. Hence the design of the estimators is done in discrete-time. Unknown Input Approach Consider the discrete-time linear stochastic system model x[k + 1] = �x[k] + �dd[k] + ![k] z[k] = Cx[k] + v[k], (7) where x is the state, d is the unknown input (disturbance) of the system and z is the measurement of the output y[k] = Cx[k]. � and �d are obtained from the matrices of Eq.(6) in the form of � = eATs , �d = Z Ts 0 eAtB dt, (8) where Ts is the sample time. The measurement and pro- cess noise v and ! are assumed to be additive, white, and Gaussian with mean of zero. Notice, that in Eq.(7), all of the system inputs are considered to be disturbances. The state and the unknown input of such a system can be estimated using an ’unknown input Kalman filter’ (UIKF) [9], if the system satisfies the conditions dim{d}  dim{y} (9) rank{C} = dim{y} (10) rank{�d} = dim{d} (11) rank{C�d} = dim{d} . (12) Additive Disturbance Approach In the case of input disturbances, load prediction is of- ten applied. Let us introduce the discrete-time nonlinear model of a system in the form x[k + 1] = �(x[k], ũ[k]) + ![k] z[k] = g(x[k]) + v[k], (13) where ũ is the disturbed input of the system, and � and g are nonlinear functions. We will follow the assumption, that the disturbance is additive to the input, thus ũ[k] = u[k] + ⇤d[k], where u is the vector of actuating signals, and ⇤ is a constant matrix. We also assume a constant disturbance model with process noise !d, d[k + 1] = d[k] + !d[k], resulting in the model x̃[k + 1] = ˜�(x̃[k], u[k]) + !̃[k] z[k] = g̃(x̃[k]) + v[k]. (14) In Eq.(14) x̃T = � xT dT � is the extended state and !̃T = � !T !T d � is the extended process noise. The ex- tended nonlinear mappings of the system are given by g̃(x̃[k]) = g(x[k]) and ˜�(x̃[k], u[k]) = ✓ �(x[k], u[k] + ⇤d[k]) d[k] ◆ . (15) In most cases the continuous-time differential equa- tion of the observed system is available, in the form 73 M uF uC r ` ✓ J ⇢ cartwinch m load rail Figure 1: The two-dimensional ovearhead crane of Eq.(4). The discrete-time state-transition equation can then be obtained by �(x[k], u[k]) = x[k] + (k+1)TsZ kTs f(x, u[k]) dt, (16) where the input is piecewise constant and u[k] = u(kTs). If the integral in Eq.(16) cannot be given in closed form and Ts is sufficiently small, the left rectangle rule approx- imation (k+1)TsZ kTs f(x, u[k]) dt ⇡ Tsf(x[k], u[k]). (17) can be used. Substituting this into Eq.(15), the expression reads ˜�(x̃[k], u[k]) = ✓ x[k] + Tsf(x[k], u[k] + ⇤d[k]) d[k] ◆ . (18) For the observation of the system described by Eq.(14) and Eq.(18) the ’unscented Kalman fil- ter’ (UKF) [10] is applicable. We chose one of the implementation variants of this filter from Ref.[10]. Application of the Methods to Overhead Cranes System Dynamics We will now consider the state estimation problem of the two-dimensional overhead crane (Fig.1). Let us denote the horizontal displacement of the cart by `, the length of the rope by r and the angle between the vertical and the rope by ✓. The generalized coordinates of the system are chosen to be q = � ` r ✓ �T , and the state is x = � ` r ✓ ˙` ṙ ˙✓ �T . Based on the sensors available, the output is y = � ` r ˙` ṙ �T . Overhead cranes is a typical example of underactu- ated systems. It is actuated by two motors, one apply- ing the force uF on the cart, and another delivering the estimated system s ⇤ + + u u system model s = ⌘(x, ẋ, u) ẋ = f(x, ũ) y = g(x) ũ y R ẋ x Figure 2: Dynamics of the estimated system torque uC to the winch. The system input consists of these, u = � uF uC �T . The parameters of the system are the mass of the cart M , the mass of the load m, the moment of inertia of the winch J and the radius of the winch ⇢. The acceleration in the gravitational field is denoted by g. The following assumptions are made of the system. • The rope connecting the load to the winch is mass- less and behaves as a rigid rod during the motion. • The load is a point mass. • M , m, J and ⇢ are known. • In the initial state of the system the unmeasured states are zero: ✓ = 0 rad and ˙✓ = 0 rad/s. • The effect of the aerodynamic resistance is negligi- ble. The Lagrangian of overhead cranes reads L = 1 2 (M +m) ˙`2 + 1 2 J ⇢2 ṙ2 ++ 1 2 m ⇣ ṙ2 + r2 ˙✓2 ⌘ +m ˙lṙ sin ✓ +m ˙lr ˙✓ cos ✓ +mgr cos ✓. (19) Using Eq.(1), the model can be written in the form of Eq.(2), where the expressions become H(q) = 0 B@ � M +m sin 2 ✓ � m sin ✓ 0 m sin ✓ ⇣ J ⇢2 +m ⌘ 0 cos ✓ 0 r 1 CA , (20) h(q, q̇) = � 0 @ m(r ˙✓2 + g cos ✓) sin ✓ mr ˙✓2 +mg cos ✓ �2ṙ ˙✓ � g sin ✓ 1 A , (21) hs(s) = � 0 @ sF sC 0 1 A , ⌧ = 0 @ uF �uC ⇢ 0 1 A , (22) where sF is the frictional force between the cart and the rail and sC is the frictional force applied to the winch. The continuous-time linear model of the crane is com- puted around the setpoint x0 = � 0 r0 0 0 0 0 �T , u0 = � 0 mg⇢ �T , and sF = sC = 0 is substituted 74 Table 1: Numeric values of the parameters Value Unit M 5 kg m 0.05 kg J 3.802 · 10�4 kg m 2 ⇢ 0.02 m g 9.81 m s �2 Ts 1 ms r0 0.47 m `D 0.5 m into Eq.(5). The resulting Jacobians are A = 0 BBBBBB@ 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 m M g 0 0 0 0 0 0 0 0 0 0 0 �m+M Mr0 g 0 0 0 1 CCCCCCA , B = 0 BBBBBB@ 0 0 0 0 0 0 1 M 0 0 � ⇢ m⇢2+J � 1 Mr0 0 1 CCCCCCA , C = 0 BB@ 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 CCA . (23) Substituting Eq.(22), the disturbance torque term in Eq.(3) reads ⌧̃ = 0 @ uF + sF � 1 ⇢ (uC � ⇢sC) 0 1 A . (24) Eq.(24) shows, that the disturbed input of the system can be introduced in the form ũ = u+ ⇤s, where ⇤ = ✓ 1 0 0 �⇢ ◆ . (25) The frictional forces usually depend on the actuating signals, state variables, and the derivative of the state as well. The general characteristics of the friction reads as s = ⌘(x, ẋ, u). The block diagram of the dynamics of the system is illustrated in Fig.2. Only the system model block is used for estimator design, which is Eq.(5) if the effect of the friction is disregarded, hs(s) = 0. Proposition 1. The discrete-time linear crane model ob- tained from the matrices in Eq.(23) and using Eq.(8) sat- isfies the design conditions of the UIKF in Eq.(9-12). Proof. It is clear that dim{y} = 4, dim{d} = 2 and rank{C} = 4, thus Eq.(9) and (10) are satisfied. Us- ing Eq.(8) and the numerical values in Table 1, �d = 10 �2 0 BBBBBB@ 0.0001 0 0 �0.025 �0.0002 0 0.2 0 0 �49.975 �0.4255 0 1 CCCCCCA , C�d = 10 �2 0 BB@ 0.0001 0 0 �0.025 0.2 0 0 �49.975 1 CCA . (26) We obtain rank{�d} = 2 and rank{C�d} = 2, which are both equal to the number of unknown inputs. This sat- isfies Eqs.(11) and (12). The continuity argument holds for the rank of the matrices in Eq.(26), therefore there ex- ists a neighbourhood around our parameter set where the design conditions are still satisfied. Simulation The estimation results were obtained by simulation of overhead cranes in a closed-loop scenario using a simple discrete-time pole-placement-based linear controller. The simulation also took into account friction. There are var- ious models of the friction phenomenon, ours included Coulomb, Stribeck, viscous friction, and stiction [11]. The numerical values of the parameters used in the simu- lation are given in Table 1. Because of friction, but mainly stiction, the position of the crane can not be accurately controlled without in- tegrators in the controller. Consequently we included the integral of the measurable positions in the controller de- sign as proposed on page 309 of Ref. [8]. The outputs to be integrated are yI = � ` r �T . In- tegrals of these quantities are approximated by the left rectangle rule as xI [k+1] = xI [k] +TsyI [k] = xI [k] +TsCIx[k], (27) where CI is defined so that yI [k] = CIx[k]. In our case CI is the first two rows of C in Eq.(23). Using Eq.(27) and Eq.(7) a new system model is introduced in the form x̃[k + 1] = ˜ �x̃[k] + ˜ �u[k] y[k] = ˜Cx̃[k]. (28) where the expanded state is x̃T = � xT xT I � and the matrices are ˜ � = ✓ � 0 TsCI I ◆ , ˜ � = ✓ �d 0 ◆ , and ˜C = � C 0 � . (29) For the controller design the process and measurement noises are omitted. 75 R DAC system in Fig. 2 ADC UKF K KI u[k] y y[k]u u0 + +- - + yID - + xD ADC x[k] x yI [k] = (`[k] r[k])T y[k] yI [k] x̂[k] ŝ[k] UIKF x̂[k] ˆ̃ u[k] u[k] CKF x̂[k] Figure 3: Closed-loop estimation setup used in the simulation 0 1 2 3 4 5 6 7 time [s] �10 0 10 20 30 40 50 60 u F , ˆ̃ u F , [N ] uF UIKF UKF 0 1 2 3 4 5 6 7 time [s] �1.0 �0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 u C , ˆ̃ u C , [N m ] uC UIKF UKF Figure 4: Actuating signals and the estimated inputs of the system We want the state-feedback controller to output u[k] = � ˜Kx̃[k] = �Kx[k] � KIxI [k]. The gain ˜K can be computed using the eigenvalue placement method implemented in Matlab’s place algorithm: ˜K = place(˜�, ˜�, �), where � contains the chosen eigenvalues of the closed-loop system. In continuous time they are �1 + 0.5j, �1� 0.5j, �5.1, �5.2, �5.3 �5.4, �5.5, and �5.6. The controller gains are K = 0 BBBBBB@ 100.594 �0.1273 0.6011 �1.6765 �185.407 0.0458 81.4011 �0.0558 0.0562 �0.3166 �4.2257 �0.0195 1 CCCCCCA T , and KI = ✓ 46.641 1.6026 �0.074 �2.955 ◆ . (30) The closed-loop setup is illustrated in Fig.3. Here the integration is done using the left rectangle rule, and the 0 1 2 3 4 5 6 7 time [s] �0.1 0.0 0.1 0.2 0.3 0.4 0.5 m ea su re d st at es [m , m /s ] � r �̇ ṙ 0 1 2 3 4 5 6 7 time [s] �0.4 �0.3 �0.2 �0.1 0.0 0.1 0.2 u n m ea su re d st at es [r ad , ra d /s ] � �̇ Figure 5: State evolution of the controlled crane using the true states as feedback ADC and DAC blocks are analogue-digital and digital- analogue converters respectively. The DAC uses the zero- order hold method for signal reconstruction. The desired state of the system is xD = � `D r0 0 0 0 0 �T , `D is the desired position of the cart and yID = CIxD. Since we only want to demonstrate the results of the state estimation, first the real states are fed back to the controller. To compare the accuracy of our me- thods to a traditional solution, the classical Kalman filter (CKF) [12] is also implemented. As a result, the Kalman filters compute the estimated state of the system x̂[k], and they also provide information about the input. The UIKF estimates the friction loaded actuating signals directly, while the UKF estimates the frictional forces superposed to the actuating signals. The actuating signals along with the estimates of the real inputs pro- vided by the filters can be seen in Fig.4. Here, the esti- mated input of the UKF is computed as ˆũ = u[k]+⇤ŝ[k]. Because of friction, the real inputs are estimated to be 76 0 1 2 3 4 5 6 7 time [s] 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 0.00035 ` � ˆ̀ [m ] UIKF UKF 0 1 2 3 4 5 6 7 time [s] �0.001 0.000 0.001 0.002 0.003 0.004 0.005 ˙̀ � ˆ̇̀ [m /s ] UIKF UKF Figure 6: Error of the position and velocity estimations much less than the actuating signals. The large offset in Fig.4 is a consequence of the stiction, the cart will not move until uF becomes greater than a certain value (in our case about 28 N). When the cart finally stops mov- ing, the fluctuation of uF tries to counter the load sway, during which the measured states do not change. This is the reason why the estimators find the input to be zero. Fig.5 shows the transient of the states in a closed- loop. Results show that the controller can bring the sys- tem into the desired state, but friction prevents it from stopping the low amplitude swinging of the load, because it makes quick, short motions impossible. Error estimates are illustrated by Figs. 6 and 7. In all cases, the error be- comes zero when the crane reaches the desired state. This is not true for the angle and its derivative in Fig.8. The error decreases at the end of the state transient, but does not become zero. When the cart stops and the load swings, the measurements are all constant thus they carry no information regarding ✓ and ˙✓. In these situations the cart and winch are usually also motionless due to fric- tion. Consequently, the estimation error of these quanti- ties never becomes zero. It is possible to decrease the estimation error using two laser slot sensors. With the help of such devices, ✓ can be measured accurately in two positions and through the application of proper sensor fusion techniques the es- timation of the angle can become more precise in between the chosen positions as well. In Fig.9 the estimation error of ✓ and ˙✓ is illustrated using CKF. The filter was designed for the dynamics given by Eq.(7) except that it does not consider the inputs unknown. The results show that the error is substantially higher than in the case of the UIKF or UKF. 0 1 2 3 4 5 6 7 time [s] �0.00005 0.00000 0.00005 0.00010 0.00015 r � r̂ [m ] UIKF UKF 0 1 2 3 4 5 6 7 time [s] �0.0010 �0.0005 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 ṙ � ˆ̇ r [m /s ] UIKF UKF Figure 7: Error of the rope length and hoisting velocity estimations 0 1 2 3 4 5 6 7 time [s] �0.003 �0.002 �0.001 0.000 0.001 0.002 0.003 0.004 0.005 � � �̂ [r ad ] UIKF UKF 0 1 2 3 4 5 6 7 time [s] �0.020 �0.015 �0.010 �0.005 0.000 0.005 0.010 0.015 0.020 0.025 �̇ � ˆ̇ � [r ad /s ] UIKF UKF Figure 8: Error of the angle and angular velocity estimations Fig.8 shows the estimation error of the UKF to be less. We applied this estimator in a closed-loop scenario identical to the one presented in Fig.3, but instead of x[k] the x̂[k] of the UKF was used in the feedback loop. The state evolution of the crane in this simulation is illustrated in Fig.10. The system’s behaviour is only slightly differ- ent from the case when the actual states were used in the feedback control in Fig.5. 77 0 1 2 3 4 5 6 7 time [s] �0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 � � �̂ [r ad ] 0 1 2 3 4 5 6 7 time [s] �5 �4 �3 �2 �1 0 1 2 3 �̇ � ˆ̇ � [r ad /s ] Figure 9: Error of the angle and angular velocity estimations using CKF Conclusion It was shown that in most mechanical systems the dis- turbing effect of friction can be reduced to the actuating signals. Based on this statement, two estimation meth- ods were provided for the computation of the unmeasured states. One of these techniques considers the inputs com- pletely unknown and uses UIKF for the estimation. The other algorithm treats the inputs as distorted by an addi- tive disturbance, and applies UKF with a load predictor extension. Simulation results in a closed-loop controlled scenario were provided to prove the applicability of the concepts. It was also pointed out, that the precision of the estimators could be improved using laser slot sensors. This will be the subject of a forthcoming paper. REFERENCES [1] HYLA P.: The crane control systems: A survey, Proc. 17th Int. Conf. Methods and Models in Au- tomation and Robotics (MMAR), 2012, 505–509 [2] SCHINDELE D., MENN I., ASCHEMANN H.: Non- linear optimal control of an overhead travelling crane, Proc. 18th IEEE Int. Conf. Control Applica- tions, 2009, 1045–1050 [3] BOUSTANY F., D’ANDREA NOVEL B.: Adap- tive control of an overhead crane using dynamic feedback linearization and estimation design, Proc. IEEE Int. Con. Robotics and Automation, 1992, 1963–1968 0 1 2 3 4 5 6 7 time [s] �0.1 0.0 0.1 0.2 0.3 0.4 0.5 m ea su re d st at es [m , m /s ] � r �̇ ṙ 0 1 2 3 4 5 6 7 time [s] �0.4 �0.3 �0.2 �0.1 0.0 0.1 0.2 u n m ea su re d st at es [r ad , ra d /s ] � �̇ Figure 10: State evolution of the controlled crane using the UKF’s estimated states as feedback [4] NEUPERT J., HILDEBRANDT A., SAWODNY O., SCHNEIDER K.: Trajectory tracking for boom cranes using a flatness based approach, Proc. Int. Joint Conf. SICE-ICASE, 2006, 1812–1816 [5] RÓZSA T., KISS B.: Tracking control for tow- dimensional overhead crane, Proc. 8th Int. Conf. Informatics in Control, Automation, and Robotics, ICINCO, 2011, 1, 427–432 [6] KISS B., LEVINE J., MULLHAUPT P.: A simple output feedback PD controller for nonlinear cranes, Proc. 39th IEEE Conference on Decision and Con- trol, 2000, 5, 5097–5101 [7] KISS B., LÉVINE J., MÜLLHAUPT P.: Modelling, flatness and simulation of a class of cranes, Electri- cal Engineering, 1999, 43, 215–225 [8] LANTOS B.: Theory and design of control systems I., Akadémia Kiadó, 2001 (in Hungarian) [9] DAROUACH M., ZASADZINSKI M., ONANA A.B., NOWAKOWSKI S.: Kalman filtering with unknown inputs via optimal state estimation of singular sys- tems, Int. J. Sys. Sci., 1995, 26(10), 2015–2028 [10] HAYKIN S.S.: Kalman filtering and neural net- works, John Wiley & Sons, Inc., 2001 [11] ARMSTRONG-HÉLOUVRY B., DUPONT P., DE WIT C.C.: A survey of models, analysis tools and compensation methods for the control of machines with friction, Automatica, 1994, 30, 1083–1138 [12] KALMAN R.E.: A new approach to linear filtering and prediction problems, Journal of basic Engineer- ing, 1960, 82, 35–45 Advertise upcoming meetings, conferences, workshops, new products or services; make public announcements, introduce your research laboratory on the pages of the Hungarian Journal of Industry and Chemistry Please contact us if interested! EDITORIAL OFFICE: UNIVERSITY OF PANNONIA P.O. BOX 158, VESZPRÉM H-8201 (HUNGARY) Tel.: +36 (88) 624-746, E-mail: hjic@almos.uni-pannon.hu; web: hjic.mk.uni-pannon.hu Felelős szerkesztő: Szilágyi Róbert Károly, PhD Kiadja: Pannon Egyetem, 8200 Veszprém, Egyetem u. 10. Levélcím: H-8201 Veszprém, Postafiók 158, Tel.: (88) 624-000 Felelős kiadó: a Pannon Egyetem, Mérnöki Kar dékánja HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 42(2) pp. 79–84 (2014) CHALLENGES OF ENTERPRISE POLICY COMPLIANCE WITH SMARTPHONE ENABLEMENT OR AN ALTERNATIVE SOLUTION BASED ON BEHAVIOUR-BASED USER IDENTIFICATION SÁNDOR DOBOS! AND ATTILA KOVÁCS Department of Information Technology, ELTE University Budapest, Budapest, 1117, HUNGARY !Email: sandor.dobos@hu.ibm.com Current trends show the intense growth in the role and importance of mobile technology (smartphones, tablets, etc.) in business due to economic, social and technological reasons. The social element drives a powerful convenience expectation called “Bring Your Own Device” (BYOD) for taking notes and accessing internal and external network resources. Apparently, the future is leading toward a more extensive enablement of smartphones and tablets with their enterprise applications. Internal security standards along with applicable regulatory ones to achieve ‘policy enforcement’ as types of solutions and controls; however, this allows for merely one aspect of compliance. An alternative solution could be behaviour-based analysis to identify the user, attacker or even a malicious program accessing resources on phone or internal networks. Complex networks can be defined by graphs, such as connections to resources on smartphones and serve as a blueprint. In case the motif is different from the user’s actual behaviour, the company can initiate specific actions to avoid potential security violations. This document reviews the IT security challenges related to smartphones as well as the concept of graph-based user identification. The challenges of the latter are the identification of motif, selection of search algorithm and defining rules for what is considered a good or bad behaviour. Keywords: mobile device security, “Bring Your Own Device” (BYOD) management, secure data communication, behaviour-based identification of threat, graph-based user identification Introduction Enterprise policy enforcement with current mobile technology management tools (Mobile Device Management (MDM), virus detection, and other necessary modules) ensure only the compliance of piece of equipment. Within enterprise, MDM Servers initiate the compliant actions, if the smartphone client fails the access privilege is taken away; however, this has an immediate impact on revenue through operational efficiency affecting the business. The enterprise compliance of the smartphone device is only a validated response, which can be altered by understanding the mobile technology, such as the application structures or the way MDM technology works. Hackers are heavily attacking smart devices with malicious software (Malware). These can be viruses, spyware, adware, and other types of attacks. A specific example is the “Obad.a Trojan”, which is now being distributed via mobile botnets. The trojans are occupying a larger space and becoming more complex, which shows the need for IT security to find new ways to detect them. The secure enterprise environment is crucial for organizations to ensure the business strategy and continuity of operation, irrespective of its environment being production, service delivery, or customer support related. Companies are following the trend to ensure efficiency and simplify service accessibility for their customers, business partners or employees. It is an essential element of the IT security strategy to be aligned with the business strategy due to (i) the intense growth of mobile technology and (ii) pressure from enterprise that has transformed business operations. Forrester Research from the Q2 Foresights Security Survey shows (Fig.1) that mobile security is of primary concern amongst enterprise leaders. CEOs are concerned about the risk of data loss, particularly due to device loss or theft. Another worry regards data protection or data leak prevention usually in connection with data related to finance and innovation. Figure 1: Degree of concern of CEOs from Forrester Research, Foresights Security Survey 80 Current Challenges Within today’s economic environment, profit-driven organizations are challenged from revenue, growth or profit target perspectives. In order to meet such expectations, organizations are under high pressure to offer innovative products (services) within newly identified channels and more accessible yet-to-be- penetrated markets to a demanding customer base, especially in the area of providing smartphone solutions. On the other hand, regulatory bodies are alarmed at mobile enablement and organizations in general striving to protect critical financial systems better, in particular via identity and access management. At the same time, organizations need to maintain operational resiliency and perform risk management assessments as an integrated part of day-to-day operations. The cost perspective is an important factor for customers facing the task of achieving a lower total cost of ownership (TCO) although, by this shift in paradigm, they move from reactive protection to proactive value creation mode. The alignment of the current IT strategy to business strategy is a challenge for corporations as IT is a service provider within an enterprise also offering business enablement leads. Companies need to manage security policy at the corporate level, assess the security health of the heterogeneous IT environment, monitor the security weakness(es) in processes and systems and, of course, meet the plethora of banking laws, regulations and standards including SOX, BASEL III, GLBA, and PCI compliance. With the adaptation of the BYOD concept in business practices, mobile technology itself contributes to additional and more complex challenges. Enterprises are not the owners of those smart devices with access to enterprise networks and their data. Fig.2 illustrates how smartphones with applications developed by unknown sources could be a threat to enterprise’s sensitive data. Business Drivers for Change As introduced above an active need is emerging for more effective mobile solutions. It is not only a reflection of the compliance requirements but a real protective solution against mobile threats without any restrictions on financial growth or cutbacks on efficiency initiatives. There are compelling reasons to act and change now including data potentially at risk, pressure from regulatory bodies and audit firms to be compliant with the applicable data handling regulations, and the implementation of mobile security into an existing security structure. Highlights of a Political-Economic-Social-Technology (PEST) Analysis Politicians are setting high standards for data management regardless of platforms, even if they are related to mobile devices or mobile applications. These regulations mostly focus on the handling of financial data such as The Sarbanes–Oxley Act of 2002 (SOX). There are new developments regarding data (pictures, text, etc.) ownership such as The Stop Online Privacy Act (SOPA). The Preventing Real Online Threats to Economic Creativity and Theft of Intellectual Property Act or in short PROTECT IP Act (PIPA) will have an even more rigorous norm for enterprises who need to meet these regulations in order to obtain permission to operate in a certain field or market (market player). The economic situation puts more pressure on businesses as the trend still shows stagnancy on a macro level while in some segments it is still in decline. Firms have high demands for business achievements measured by revenue and profit perspectives to meet shareholder expectations. One way to do so is to improve cost structures with mobile enablement (particularly BYOD). Society wants to live more comfortably by searching for easily accessible services, such as online mobile- enabled bank transfers or online mobile retail orders. Companies aim to reach their customers via mobile applications that are more easily accessible through user-friendly interfaces, such as social network platforms. As usual, this trend is typical for major market areas, such as North America or, Western Europe; however, other regions are expected to follow. The trend in technology continues to be more vivid from a computing power and portability perspective. Internet Service Providers are granting access to higher bandwidth networks for new applications and online content. Policy Enforcement and its Weaknesses To manage the various types of mobile devices in an enterprise, a solution was introduced in 2008 by STRICKLEN, MCHALE, CAMINETSKY and REDDY based on Pattern Definition for Mobile Device Management (MDM). MDM is the most common approach how mobile devices (including smartphones, tablets) are managed within the enterprise as a management tool for monitoring and policy enforcement. Even though this technology is available and wide-spread, according to ABI research, mobile threats grew by 261% just in the Figure 2: Graphical illustration of mobile phone technology security threats 81 last quarter of 2012; however, MDM technology is not growing at a comparable rate to the number of vulnerabilities and attacks. One of the components of MDM is the registration module, which is responsible for identifying mobile devices within the enterprise network, where the identifier indicates the platform of each mobile device, e.g. smartphone and tablet. The other significant element is the management module. Its job is to receive a management function definition performed on mobile devices using identifiers from mobile devices. The management module is responsible for instructions for the first platform to perform the management tasks and provides instructions for the other modules specific to the second platform to perform the management task on at least the second device. However, MDM is not the ultimate solution for mobile device security within the enterprise environment due to some concerns regarding the verifiable device integrity. There are attempts at identifying MDM solutions to detect modifications of the underlying platform, but since the MDM agent has limited privileges and was susceptible to compromise by malicious privileged software, these stand little chance of detecting a targeted attack. The US National Security Agency (NSA) describes an immutable cryptography, as the ‘root of trust’ on a specific platform, to be available for leveraging by MDM or other software, providing a means of countering this type of threats. Encryption enables devices to attest their integrity on an enterprise and carry out any local policy decisions. The availability of this ‘root of trust’ to other software can vigorously complement a chain of trust that begins when booting and extends into the system runtime. An additional benefit of an immutable root of trust is that it allows a company to bind the unique identifier of that device together with other credentials to restrict company access to only those devices. In effect, the device itself can become one of the factors of a multi-factor access. In summary, MDM’s current capabilities support the BYOD and enterprise-owned use of mobile devices with certain gaps for high-security issues. Management capabilities are limited to those provided for MDM products by the underlying mobile OS (iOS, Android, and other platforms), and therefore, these capability gaps cannot be closed by MDM providers alone. On- going cooperation between enterprise customers, OS vendors, and MDM vendors is critical to the continued advancement of enterprise-level security for mobile devices. Closing stated holes will enable the deployment of commercially available mobile devices to tackle high-security use issues common in sensitive industrial and governmental environments. Secure data networks need to guarantee integrity, confidentiality and availability; this is when security is fulfilled. Policy enforcement with the MDM can be too liberal, but then there is no real need for MDM or strictly constraining users thus limiting the value-add. The proposed solution is to identify the user based on behaviour and compare behaviour changes. Behaviour Analysis and Challenges Growing recognition drove the importance of network science related to the behaviour of sophisticated systems that is shaped by relations among their constituent elements. The rising availability and tractability of large and high quality data sets on a wide range of complex systems [1-3] have led to a primary insight: substantially diverse complex systems often share core key organizational principles. These can be quantitatively characterized by the same parameters, which means that they show remarkably similar macroscopic behaviour despite reflective differences in the low-level details of the components of each system or their mechanisms of interaction. The behaviour as described above can be modelled by mathematical graphs, where the graphs are defined simply as a set of nodes called vertices linked together in direct or indirect ways by connections (edges). From a mathematical point of view, G(V, E) are canonical graphs as the vertices and edges are labelled. In network science, methodological advances permit research to quantify other topological properties of complex systems, such as modularity [4], hierarchy [5], centrality [6] and the distribution of network hubs [7,8]. There have also been significant efforts to form the development or evolution of complex networks [9], to link network topology to network dynamics, and to explore network robustness and vulnerability. These topics are likely to become more relevant in relation to behavioural studies. Structural and functional behaviour maps can be created using graph theory through the following four steps: 1. Define the network nodes. 2. Estimate a continuous measure of association between nodes. 3. Generate an association model by compiling all pairs of associations between nodes and usually apply a threshold to each element of this model to produce a binary adjacency matrix or undirected graph. 4. Calculate the network parameters of interest in this graphical model of the behaviour network and evaluate them against the equivalent parameters of a population of random networks. The elements such as defined network nodes, relationship of nodes, generated association matrix and network parameters give the basis of the behaviour type of “blueprint” of a user. A mobile user can be defined as an identity based on certain features of this behavioural “blueprint” as the key points of a fingerprint within the enterprise company. In Fig.3 the darker highlighted edges could be the essential elements for the identifier. Access to resource driven networks is defined in the graph theory as a set of nodes or vertices and the edges or lines, plus the connections between them. The graph topology can be quantitatively captured by a wide variety of measures used to generate the key points of the blueprint. The most important measurement is the node degree. The degree of a node is the number of 82 connections it links to the rest of the network. Based on this definition, it is the most elementary network measure, and most other measures are ultimately linked to the node degree. The degrees of all the network’s nodes identified form a degree distribution [10]. In randomly selected networks all of the connections are equally probable, i.e. the result is a Gaussian distribution. However, complex networks in general have non-Gaussian distributions, often with a long tail towards higher degrees. The degree distributions of the scale-free networks follow a power law [11]. This landmark study [11] was the first to describe the scale- free organization of many complex networks and proposed a simple growth rule for their formation. The number of connections could give a large graph network that could be used for usability needs to define clustering, coefficients, and motifs. If the nearest neighbours of a node are also directly linked to each other, they form a cluster, and the clustering coefficient quantifies the number of connections that subsist between the nearest neighbours of a node as a quantity of the maximum number of possible connections [12]. The randomly selected networks exhibit low average clustering, whereas compound networks displays high clustering. Interactions between neighbouring nodes can also be counted by counting the occurrence of small motifs of interconnected nodes [13]. The distribution of different motif classes (nodes EFL and KMJ in Fig.3) in a network provides information about the types of local interactions that the network can support [14]. The numbers of reachable nodes and hops needed for connectivity are defined by the path length and efficiency. The path length is the minimum number of edges that must be passed through to start moving from one node to another. Random, complex networks have short mean path lengths; they exhibit high global efficiency of parallel information transfer, whereas regular patterns exhibit long mean path lengths. Efficiency is inversely related to path length. Nevertheless it is numerically easier to estimate topological distances between elements of disconnected graphs. Moreover, the link density or cost of the route provides further description. The connection density is the actual number of the network edges in the graph as a proportion of the total number of possible edges and is the simplest estimator of the physical cost of the network. Hubs are nodes of high degree, where the centrality of a node measures how many of the shortest paths exist between all other node pairs in the network passing through it. It can be asserted that nodes with high centrality are thus crucial to efficient communication [15]. The importance of an individual node to network effectiveness can be assessed by deleting it and estimating the efficiency of the severed network. The robustness property refers either to the structural integrity of the network following the deletion of nodes, edges or effects of perturbations on local or global network states. Finally, modularity is the property that has significant influence on the “blueprint” of the behaviour. There are algorithms that estimate the modularity of the network on the basis of hierarchical clustering [16]. Each module contains several densely interconnected nodes, and there are relatively few connections between these nodes in different modules. Airline hubs described this as a function of their roles in community structure [17]. Provincial hubs are connected mainly to nodes in their modules where the connector hubs are connected to nodes in other modules. Graph Based Analysis of Smartphone Usage Previously, the challenges of using mobile devices in the enterprise environment were discussed. It is important to identify the activities on the device in order to take either preventive or corrective actions. G(V, E) used for the behaviour identification and description graphs. Vertices are defined as applications, more precisely, vertices within the device or network addresses and sites. Edges are the connections between the vertices applications called websites, network addresses, user- or process-initiated connections. • Network nodes can be defined as the resource access, i.e. applications, programs, and device elements access initiated by the user or application calls. • Association between network nodes could be defined as the connections of the applications (calls, connections, as web access), usually resources are connecting together. • Based on the above elements, network nodes and associations provide the details to calculate the association matrix. • The following steps include network parameter analysis calculated from the above. The question is not whether mobile activity can be defined as a complex network and describable by graph, but whether the process identies a user by a motif. However, before this question can be answered certain challenges need to be faced: • identify motifs as maximum independent set, • search motif in the complex network, • define „approved behaviours” according to the identification of user behaviour. Figure 3: An example of a mobile user’s identity on the basis of a “blueprint” network 83 Identification of the Motif in the Graph The task is to find the Maximum Independent Set (MIS) of a particular graph as was defined by GAREY and JOHNSON [18] as NP-complete and remains so even for bounded degree graphs. According to FEIGE et al. [19] MIS cannot be approximated even within a factor of |V|1 − o(1) in polynomial time. The greedy algorithm (GMIS) can provide a solution for the identification of the maximum independence set. The algorithm selects a vertex of minimum degree, deletes that vertex and all of its neighbours from the graph, and repeats this process until the graph becomes empty. A recent detailed analysis of the GMIS algorithm has shown that it produces reasonably good approximations of the MIS for bounded- and low- degree graphs defined by HALLDORSSON and RADHAKRISHNAN [20]. In particular, for a graph G with a maximum degree Δ and an average degree !, the size |I| of the MIS satisfies Eq.(1). |!| ≤ min! ∆!! ! |GMIS ! |, !!!! |GMIS ! | (1) where |GMIS(G)| is the size of the approximate MIS found by the GMIS algorithm. Eq.(1) provides an upper-bound of the number of edge-disjoint embeddings of a particular sub-graph, and will use this bound to obtain a computationally tractable problem formulation that is guaranteed not to miss any sub-graphs that can potentially be frequent. Searching for a Discover Motif in the Graph The GROCHOW-KELLIS method [21] can be used to find all examples of sub-graphs of a given size, similar to exhaustive methods. In the background, all non- isomorphic graphs are generated of a particular size using MCKAY’S tools [22]. Then for each graph, the method evaluates its significance. Due to symmetries, a set sub-graph of G may be mapped to a set query graph H multiple times. Therefore, a simple mapping-based search for a query graph will locate every case of the query graph as many times as the graph has symmetries. To avoid this, the method computes and enforces several symmetry-breaking conditions, which ensure that there is an exclusive map from the query graph H for each case of H in G, so that the search only spends time finding each instance once IsomorphicExtensions (f,H,G[,C(h)]: Finds all isomorphic extensions of partial map f : H → G [satisfying C(h)] Start with an empty list of isomorphism. Let D be the domain of f. If D=H, return a list consisting solely of f. (Or write to disk.) Let m be the most constrained neighbour of any d � D (constrained by degree, neighbours mapped) For each neighbour n of f (D) If there is a neighbour d � D of m such that n is not neighbours with f(d), or if there is a non-neighbour d � D of m such that n is neighbours with f(d) [or if assigning f (m)=n would violate a symmetry- breaking condition in C (h)], then continue with the next n. Otherwise, let f=f on D, and f (m)=n. Find all isomorphic extensions of f. Append these maps to the list of isomorphism. Return the list of isomorphism The symmetries of a graph H equals automorphism (self-isomorphism). The collection of automorphism of H is indicated by Aut(H). For a set A of automorphism, two nodes are stated to be “A-equivalent” if there is some automorphism in A which maps one to the other, or just “equivalent” if A =Aut(H). Given a set of conditions C, α preserves the conditions C if, given a labelling L1 of H which satisfies C, the corresponding labelling L2: H → Z given by L2(n) = L1(α(n)) also satisfies C. We are thus searching for setting C, such that the only automorphism preserving C is the identity. This ensures exactly one map from H onto each of its instances in G to satisfy the conditions. To find these conditions, an Aut(H)-equivalence class {n0, ..., nk} of nodes of H, and the condition L(n0) < min (L(n1), ..., L(nk)) imposed. Any automorphism must send n0 to one of the ni, since these are all of the nodes equivalent to n0. However to preserve this state, an automorphism must send n0 to itself. Then the process continues recursively, replacing Aut(H) with set A of automorphisms that send n0 to itself. Because FindSubgraphInstances starts with a particular node that node can be considered already fixed. The main reasons for the selection of the GROCHOW and KELLIS method are summarized as follows: • Capable of finding more significant motifs by enabling exhaustive discovery of motifs up to seven nodes. To find even larger motifs, the method samples a connected subgraph, and then finds all its occurrences and assesses their significance using this method. This practice has enabled the algorithm to find motifs of up to fifteen nodes and examine subgraphs of up to thirty-one nodes; • Capable of quering a particular subgraph by querying whether a particular subgraph is significant; • Capable exploring motif clustering; because the algorithm finds all occurrences of a given subgraph, it can be used to examine how these cases cluster together to form larger structures; • Time and space applied to all subgraphs of a set size, takes exponentially a smaller amount of time than previous methods, even when implementing the previous method with the hashing scheme. Conclusions Apparently the current policy enforcement solutions are not sufficient because they mainly focus on policy compliance based on the response of the device. The installation of malicious applications could alter the mobile devices response (jailbreak on IOS or Android 84 and software imitating a response to the MDM server) taking to more serious security issues. Similarly to the recently identified malware Oldboot.B, which can install malicious applications in the background, it can inject malicious modules into the system process that prevents malware applications from uninstalling. Oldboot.B can change the browser setup (set a new and unwanted home page), and it also can uninstall or disable installed Mobile Antivirus software, and even steal data such as credit card information or any other critical data. This malware is especially dangerous as it implements evasion techniques to stay undetected. In order to optimize current threat detection, which is also in line with business requirements, additional and different methodologies are required such as user identification based on the blueprint of the behaviour of the user. To create a working model of the graph-based user identification, several tasks await completion. The running tasks are the edges started by the system or user defined as the connections between them. Once the complex network and graph are set up the motif can be identified with a greedy algorithm and search of the global graph with the GROCHOW-KELLIS Algorithm. In order to have a behaviour-based decision making system, all of the above algorithms will be tested to further develop the behaviour-based system. REFERENCES [1] AMARAL L.A.N., SCALA A., BARTHELEMY M., STANLEY H.E.: Classes of small-world networks, Proc. Natl Acad. Sci. USA 2000, 97, 11149–11152 [2] AMARAL L.A.N., OTTINO J.M.: Complex networks. Augmenting the framework for the study of complex systems, Eur. Phys. J. B 2004, 38, 147– 162 [3] BARABÁSI A.L., OLTVAI Z.N.: Network biology: understanding the cell’s functional organization, Nature Rev. Genet. 2004, 5, 101–113 [4] GIRVAN M., NEWMAN M.E.J.: Community structure in social and biological networks. Proc. Natl Acad. Sci. USA 2002, 99, 7821–7826 [5] RAVASZ E., BARABÁSI A.L.: Hierarchical organization in complex networks, Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 2003, 67, 026112 [6] BARTHÉLEMY M.: Betweenness centrality in large complex networks, Eur. Phys. J. B 2004, 38, 163– 168 [7] GUIMERA R., AMARAL L.A.N.: Functional cartography of complex metabolic networks, Nature 2005, 433, 895–900 [8] GUIMERA R., MOSSA S., TURTSCHI A., AMARAL L.A.: The worldwide air transportation network: anomalous centrality, community structure and cities’ global roles, Proc. Natl Acad. Sci. USA 2005, 102, 7794–7799 [9] KASHTAN N., ALON U.: Spontaneous evolution of modularity and network motifs, Proc. Natl Acad. Sci. USA 2005, 102, 13773–13778 [10] AMARAL L.A.N., SCALA A., BARTHELEMY M., STANLEY H.E.: Classes of small-world networks, Proc. Natl Acad. Sci. USA 2000, 97, 11149–11152 [11] BARABÁSI A.L., ALBERT R.: Emergence of scaling in random networks, Science 1999, 286, 509–512 [12] WATTS D.J., STROGATZ S.H.: Collective dynamics of “small-world” networks, Nature 1998, 393, 440–442 [13] AMARAL L.A.N., OTTINO J.M.: Complex networks. Augmenting the framework for the study of complex systems, Eur. Phys. J. B 2004, 38, 147– 162 [14] BARABÁSI A.L., OLTVAI Z.N.: Network biology: understanding the cell’s functional organization. Nature Rev. Genet. 2004, 5, 101–113 [15] FREEMAN L.C.: A set of measures of centrality based on betweenness, Sociometry 1977, 40, 35– 41 [16] GIRVAN M., NEWMAN M.E.J.: Community structure in social and biological networks. Proc. Natl Acad. Sci. USA 2002, 99, 7821–7826 [17] GUIMERÀ R., MOSSA S., TURTSCHI A., AMARAL L.A.: The worldwide air transportation network: anomalous centrality, community structure and cities’ global roles. Proc. Natl Acad. Sci. USA 2005, 102, 7794–7799 [18] GAREY M.R., JOHNSON D.S.: Computers and Intractability: A Guide to the Theory of np- completeness. New York: W. H. Freeman and Company, 1979 [19] FEIGE U., GOLDWASSER S., LOVASZ L., SAFRA S., SZEGEDY M.: Approximating clique is almost NP complete, Proc. 32nd IEEE Symposium on Foundations of Computer Science (FOCS), 1991, 2–12 [20] HALLDORSSON M.M., RADHAKRISHNAN J.: Greed is good: Approximating independent sets in sparse and bounded-degree graphs, Algorithmica, 1997, 18(1), 145–163 [21] GROCHOW J.A., KELLIS M.: Network Motif Discovery using Subgraph Enumeration and Symmetry-Breaking, Computer Science and AI Laboratory, M.I.T. Broad Institute of M.I.T. and Harvard, 2007 [22] MCKAY B.D.: Isomorph-free exhaustive generation, J. Algorithms, 1998, 26, 306–324 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 42(2) pp. 85–89 (2014) MEASUREMENT-BASED MODELLING AND SIMULATION OF A HYDROGEN- GENERATING DRY CELL FOR COMPLEX DOMESTIC RENEWABLE ENERGY SYSTEMS ATTILA GÖLLEI,! PÉTER GÖRBE, ATTILA MAGYAR, AND LÁSZLÓ NEUKIRCHNER Department of Electrical Engineering and Information Systems, University of Pannonia, Egyetem u. 10., Veszprém, 8200, HUNGARY !E-mail: golleia@almos.uni-pannon.hu Nowadays, the growing need for energy from renewable sources and growing revulsion towards fossil and nuclear fuels puts sustainable and green energy in the limelight. Producing (electrical) energy in domestic power plants from renewable sources (mainly solar and wind) hardly results in difficulties, but the storage of energy not consumed immediately is a great engineering challenge. In the present paper a complex model has been developed by investigating renewable energy sources, the surplus energy not actually consumed and stored in electrical vehicle (EV) batteries, the conversion to hydrogen for storage purposes and how the main grid is fed. A measurement-based model of a hydrogen generating cell developed for the simulation of complex energetic systems. The parameter estimation of the static model was based on the collected measurement data coming from the detailed examination of a built demonstration cell. The novel element of this work is the Matlab Simulink model for the hydrogen generation cell. Using this model, a dynamic simulator of a complex domestic power plant is made available using renewable energy sources and hydrogen generation cells. Hydrogen generation enables the lossless long-term storage of surplus electric energy collected, but not consumed or injected into the low voltage grid. The generated hydrogen can be consumed for transportation purposes in suitable vehicles or it can be applied in fuel cells generating direct electrical energy for energy-deficient low voltage network situations. Energetic situations potentially occurring in practice were simulated in our complex model. Simulations showed that the presented model is suitable for domestic scale low voltage complex energetic systems. Keywords: hydrogen generation, renewable energy sources, domestic power plants, modelling and simulation, measurement-based modelling Introduction Producing hydrogen gas (H2) from excess energy is not a new idea. This is an alternative way to store and convert renewable energy for further utilization. The produced hydrogen can be stored or used in power cells to be converted back to electric energy or in vehicles for hydrogen propulsion [1]. Although the described procedure is efficient and able to produce a high quantity of H2 it is not suitable for application in combination with domestic power plants. The most relevant from of H2 production is when the energy consumption and quantity of produced H2 are controlled. When power consumption and generation are continuous (and not necessarily deterministic) functions of time, the H2 production depends solely on the excess energy of the grid. The best solution is the usage of Supervisory Control and data Acquisition Systems (SCADA) of management [2]. The domestic applicability of this technology in the future depends on the cost of SCADA system installation. Producing hydrogen and oxygen gases from water using electricity in a laboratory is a simple electrochemical process that can be performed easily and in a very demonstrative way. Producing hydrogen on a large scale or in industrial quantities calls for an optimized or near-optimized cell model. In an energy demanding process only a few percent of variance in efficiency could mean a significant energy surplus or shortage [3]. The electrochemical parameters of a dry cell (Fig.1) that are used here are discussed to simulate hydrogen and oxygen gas production. Compared to wet oxyhydrogen (HHO) cells where the entire unit is underwater, the plates of dry cells are separated with rubber seals. These seals stop the water from leaking from the cell. The electrical connections and edges of Figure 1: The theoretical setup of a dry HHO block 86 the plates do not touch the electrolyte. These parts of the unit stay dry, thus the name dry cell. To make sure the gas made from the electrolyte gets out of the cell and the solution flows between the plates, there are holes on the top (for the gas) and bottom (for the electrolyte) on the metal slats (Fig.1). The application of dry HHO units has two main advantages. The surface of the dry cell plates enables one to use smaller amounts of electrolyte compared to with wet cells; therefore, the volume and weight of the cell is smaller. Furthermore, the connectors of dry cells remain dry, i.e. they do not corrode unlike to wet cells, where the connectors are underwater therefore their surface slowly corrodes [4]. Electrochemical Foundations Electrochemical cells can be considered as galvanic batteries where the electrochemical reactions are supported by an external current supply. They are composed of two electrodes and a conductive electrolyte fluid. If the electrode material does not participate directly in the electrode reaction, it is called an indifferent electrode (e.g. graphite). During electrolysis, if there is more than one possible type of electrochemical reaction, then a simple anion will detach from the positive anode (e.g. chloride), without this anion, HO- will be created by water splitting. The dissolution voltage of water is 1.23 V at 25 °C, the temperature coefficient is -0.85 mV/K, which means that at 100 °C this voltage decreases to 1.17 V. Therefore, in the light of these data, the specific energy demand to make hydrogen via electrolysis at 25°C can be calculated from Eqs.(1–4). The amount of charge needed to evolve 1 kg of H2 gas is q = zFM = 2⋅96487⋅0.5 = 96487 A s mol-1 = = 26801 Ah kg-1 (1) wH2 = qEMF = 26801⋅1.23 = 32966 Wh kg-1 (2) Since the volume of 1 kg of standard state H2 is 12474 dm3, the amount of energy required to produce 1 dm3 of H2 gas is: wH2 = 32966 12474 = 2.64 Wh dm−3 (3) To generate 1 dm3 of hydrogen gas, 1.5 dm3 of HHO gas is needed and thus the energy demand of producing 1 ldm3 of HHO gas (0.667 dm3 H2) is: wH2(HHO) = 0.667⋅2.64 = 1.76 Wh dm-3 (4) The unit has been measured at 10 different electrolyte concentrations, using different currents. At the same time, the voltage on the plates and amount of gas produced by electrolysis has also been measured. The HHO Cell Unit The setup of one block of the unit is shown in Fig.2. Usually five cells make up one block giving one gas- producing block. The block’s electrical connections are on the ends of two plates (Fig.1). Four of the six electrode plates are neutral electrodes, as there is no voltage connected to them. The potential is divided between the neutral plates according to voltage division in series connections. It means that the voltage between two electrodes is one fifth of the voltage on one whole block. In the experiment, a unit with three blocks connected in parallel has been used. Besides the HHO cell, a water reserve tank to infuse the electrolyte into the cell was necessary. A tube between the gas outlet and the tank has also been installed since due to bubbling, electrolyte comes out of the tube that needs to be recycled back into the system. Then, as the electrolyte drips back into the tank, the gas can escape into the bottle through another hose. The produced H2 volume and the production speed are measured with this bottle. A power supply (Manson SPS9600) has been connected to the electrical connections of the HHO unit, in this way the input current was controlled (Table 1). Matlab Model of the Dry Cell The model of the dry cell considered was implemented in Matlab Simulink using the SimPowerSystems Toolbox. Two unknown functional relationships between the generated H2 volume, the cell current and the KOH concentration and between cell voltage, cell current and KOH concentration were approximated using fourth and third order polynomials, respectively Figure 2: The setup of a HHO gas generator cell block Table 1: Experimental results Electrolyte concentration, g dm-3 MMWa, cm3 min-1 W-1 Gas production, dm3 min-1 Power of unit, W 1 2.13 0.20 10.8 2 2.66 0.75 34.4 3 2.66 1.37 55.8 4 2.59 1.51 82.2 5 2.72 1.90 90.6 6 2.63 2.52 119.5 7 2.67 2.96 140.0 8 2.65 2.76 125.0 9 2.46 2.28 105.6 10 1.82 2.15 103.2 a millilitres per minute per watt 87 using the Matlab Surface Fitting Tool. As the fitted polynomials do not have a physical connection to the given device, the model is applicable to any similar electrochemical H2 generation device with an electric two-pole system. In the different linear and non-linear physical and chemical models different coefficients become dominant. The voltage relationship is given by Eq.(5), where icell denotes the cell current and cKOH stands for the KOH concentration. Parameters can be found in Table 2. ucell(icell,cKOH) = p00 + p10 icell + p01 cKOH + p20 icell 2 + p11 icellcKOH + p02cKOH 2 + p30icell 3 + p21icell 2 cKOH (5) The volume of the generated H2 is given by Eq.(6). H2 (icell,cKOH ) = p00 + p10 icell + p01 cKOH + p20 icell 2 + p11 icellcKOH + p02 cKOH 2 + p30 icell 3 + p21 icell 2 cKOH + p12 icellcKOH 2 + p03 cKOH 3 + p40 icell 4 + p31 icell 3 cKOH +p22 icell 2 cKOH 2 + p13 icellcKOH 3 (6) Table 3 and Figs.3-4 show representative results for the model. As expected, the H2 generation speed decreases and the cell finally stops working as the amount of water decreases and the KOH concentration increases. A Simulink block scheme of the cell model is depicted in Fig.5. This Simulink model was validated by considering a system with the same parameters as the layout of the experimental cell. In this layout, we ran a simulation for 24 h using this model, decreasing water and increasing KOH concentrations. The results of this simulation can be seen in Fig.6. It can be seen that the hydrogen gas generated is reduced because of the rising KOH concentration. The exact values are in good agreement with our measurements. Dry Cell Model in Complex Energetic Systems The model for H2 generating cells described in the previous section was investigated in the Matlab Simulink simulation environment that studies the energy flow conditions of a complex energetic system consisting of a renewable source with a grid- Table 2: Coefficients of the polynomial relationship describing the cell voltage Value Value Value p00 1.429 p10 0.2548 p01 -0.1226 p20 -0.008571 p11 -0.01191 p02 0.008257 p30 0.0001141 p21 -8.76e-05 p12 0.0009697 Table 3: Coefficients of the polynomial relationship for the generated H2 gas Value Value Value p00 -0.1695 p10 0.1687 p01 -0.01765 p20 -0.007486 p11 -0.03234 p02 0.03446 p30 -0.0001077 p21 0.00412 p12 -0.004094 p03 -0.004061 p40 -4.269e-06 p13 6.169e-05 p22 -0.0005518 p13 0.0009544 Figure 3: Simulation of the cell model with a constant current of 5 A for 1 day Figure 4: Generated H2 as a function of KOH concentration and dry cell current Figure 5: Matlab Simulink model of the HHO cell. The functional blocks implementing Eqs.(5) and (6) are denoted by different background colours Figure 6: Simulation of a cell model with a constant current of 5 A for 1 day 88 synchronized inverter, a low voltage grid, an intermediate voltage controller [3,5] and a lithium ion battery. We replaced the lithium ion battery in this cell model, which reduces the potential energy flow modes, because this cell can only adsorb current for storing energy in hydrogen production. It cannot reverse the electrochemical process for electrical energy generation from hydrogen gas. The structure of the system can be seen in Fig.7, where it is apparent that the cell model is connected directly only to the grid-synchronized inverter module of the system. The system depicted in Fig.7 operates in different discrete states according to the energy flow direction. Four cases can be defined: • Normal inverter mode: The energy flows from the renewable source to the grid only (Fig.8A). • Normal inverter and hydrogen generation mode: The energy flows from the renewable source to both the dry cell and the grid (Fig.8B). • Hydrogen generation only mode: The energy flows from the grid to the dry cell only (Fig.8C). • Distortion reduction only mode: The energy flows from the grid into the intermediate capacitance and from the intermediate capacitance into the grid. The energy balance is zero for a period, and the active power is zero (Fig.8D). Model verification was performed by changing the energy flow modes in subsequent time intervals, and this was implemented by changing the energy balance of the system with outer current loads (Iouter load). The different values for the simulations as parameters can be Figure 7: Simulink model of a complex energetic system with an HHO cell model inside Figure 8: Complex energetic system energy flow modes: A: normal inverter mode, B: inverter and hydrogen generator mode, C: hydrogen generation only mode, D: distortion reduction only mode 89 seen in Table 4. The simulation results are shown in Fig.9, where Uconn is the effective value of the voltage at the connection point, IHHO is the current value of the dry cell, VH2 is the volume of generated hydrogen gas, cKOH is the KOH concentration of the electrolyte and Vwater is the volume of water inside the cell system. These values are plotted as a function of time. The results of the simulation show that the behaviour of the simulated electronic two-pole system is identical to that of the measured database. Conclusion We developed a complex model to investigate renewable energy sources, for the conversion of surplus energy to hydrogen gas for storage and to represent how a main grid is fed. We built a measurement-based model of a hydrogen-generating cell for the simulation of complex energy systems in MATLAB SIMULINK environment. We estimated the parameters of the model based on measurements collected during the detailed examination of a demonstration cell. We carried out a series of experiments on a HHO gas producing dry cell to find the optimal electrolyte concentration, current value, etc. or change the setup by altering the distance between the plates with KOH electrolyte solution. We monitored the experimental setup in several regards, for example cell voltage, and gas production. The novel element is the temperature and concentration dependent Matlab Simulink model of the hydrogen generation cell, which was found to be suitable for simulation purposes. We tested it in a simulation of a complex domestic power plant using a renewable energy source and hydrogen generation cell. Hydrogen generation enables the long-term storage of surplus electrical energy collected, but not consumed or injected into the low voltage grid. The generated hydrogen can be consumed by vehicles for transportation purposes or it can be applied in fuel cells generating direct electrical energy for energy-deficient low voltage network situations. We simulated all the potential energetic situations in this complex model of an energetic system. The simulations showed that the presented model of a hydrogen- generating cell performed well. ACKNOWLEDGEMENT We acknowledge the financial support of this work by the Hungarian State and the European Union under the TAMOP-4.2.2.A-11/1/ KONV-2012-0072 project. REFERENCES [1] KOUTSONIKOLASA D.E., KALDISA, S.P., PANTOLEONTOSB, G.T., ZASPALISAC, V.T., SAKELLAROPOULOSC, G.P.: Techno-Economic Assessment of Polymeric, Ceramic and Metallic Membranes Integration in an Advanced IGCC Process for H2 Production and CO2 Capture, Chem. Engng. Trans. 2013, 35, 715-720 [2] ZIOGOUA C., ELMASIDESB C., PAPADOPOULOUCA S., VOUTETAKISA S.: Supervisory Control and Unattended Operation of an Off-Grid Hybrid Power Generation Station with Hydrogen Storage Chem. Engng. Trans. 2013, 35, 529-534 [3] GÖRBE P., MAGYAR A., HANGOS K.M.: Line Conditioning with Grid Synchronized Inverter's Power Injection of Renewable Sources in Nonlinear Distorted Mains, Proc. 10th International PhD Workshop on Systems and Control, 2009, 978-980 [4] AL-ROUSAN A.A.: Reduction of fuel consumption in gasoline engines by introducing HHO gas into intake manifold, Int. J. Hydrogen Energy, 2010, 35, 12930-12935 [5] GÖRBE P., MAGYAR A., HANGOS K.M.: THD Reduction with Grid Synchronized Inverter’s Power Injection of Renewable Sources, Proc. 20th Int. Symp. Power Electronics, Electrical Drives, Automation and Motion (SPEEDAM) 2010, 1381- 1386 Figure 9: Simulation results of a complex energetic system using a cell model short time range (5 sec) Advertise upcoming meetings, conferences, workshops, new products or services; make public announcements, introduce your research laboratory on the pages of the Hungarian Journal of Industry and Chemistry Please contact us if interested! EDITORIAL OFFICE: UNIVERSITY OF PANNONIA P.O. BOX 158, VESZPRÉM H-8201 (HUNGARY) Tel.: +36 (88) 624-746, E-mail: hjic@almos.uni-pannon.hu; web: hjic.mk.uni-pannon.hu Felelős szerkesztő: Szilágyi Róbert Károly, PhD Kiadja: Pannon Egyetem, 8200 Veszprém, Egyetem u. 10. Levélcím: H-8201 Veszprém, Postafiók 158, Tel.: (88) 624-000 Felelős kiadó: a Pannon Egyetem, Mérnöki Kar dékánja HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 42(2) pp. 91–95 (2014) QUASI-POLYNOMIAL REPRESENTATION-BASED CONTROL OF MECHANICAL SYSTEMS LÁSZLÓ NEUKIRCHNERB AND ATTILA MAGYAR Department of Electrical Engineering and Information Systems, University of Pannonia, Egyetem u. 10., Veszprém, 8200, HUNGARY BE-mail: neu.groszber@gmail.com A simple kinematic model of a differential steering mobile robot is investigated using a nonlinear technique based on the quasi-polynomial representation of the dynamical model. Dynamical systems can be embedded in the generalized Lotka-Volterrs (or quasi-polynomial) form under mild conditions. Quasi-polynomial systems are good candidates for a general nonlinear system representation, since their global stability analysis is equivalent to the feasibility of a linear matrix inequality. The stabilizing quasi-polynomial state feedback controller design problem is equivalent to the feasibility of a bilinear matrix inequality. The classical stabilizing state feedback problem for quasi-polynomial systems was extended with the ability of tracking time-dependent reference signals. It is shown that the stabilizing quasi-polynomial controller design is equivalent to a bilinear matrix inequality. The results are applied to the model of the differential steering mobile robot. The goal reaching quasi-polynomial controller is shown to be a special kind of proportional state feedback. Keywords: quasi-polynomial, robotics, differential drive robot, control systems, kinematics, Lotka-Volterra system Introduction Trajectory and goal tracking of mobile robots is an inten- sively studied field of modern robotics as well as mod- ern control theory. Several papers deal with an adaptive output feedback approach [1]. On the other hand some groups try to apply neural network-based methods for the task [2]. Another direction is to describe the prob- lem as an optimal control problem and apply optimal control results for it [3]. The class of quasi-polynomial (QP) systems plays an important role in the theory of nonlinear dynamical systems because nonlinear systems with smooth nonlinearities can be transformed into quasi- polynomial form Ref.[4]. This means, that any applicable method for quasi-polynomial systems can be regarded as a general technique for nonlinear systems [5]. The aim of this work is to widen the applicability of the quasi- polynomial representation-based methods to the class of mechanical systems, more precisely to mobile robots. The goal reaching problem of a differential steering mo- bile robot is reformulated as a globally stabilizing feed- back design problem in the quasi-polynomial framework. Basic Notions Differential Drive Mobile Robot Kinematics The chosen mechanical system is the kinematic model of a two wheeled differential drive mobile robot Eq.(1). The model deals with the geometric relationships that govern the system. It calculates the motion without considering the affecting forces. The system’s states are the Cartesian coordinates x, and y and the orientation ✓ of the mobile robot. The basic kinematic model of the differential drive robot is given by ˙✓ = r 2a (!left � !right) ẋ = r 2 cos(✓)(!left + !right) ẏ = r 2 sin(✓)(!left + !right), (1) where a is half the shaft’s diameter, r is the radius of the wheels and ! is the angular velocity of the right or left wheel. To order the robot to reach a specific goal, it is acceptable to design a proportional controller to govern the expected trajectory Eq.(2) [6] . In this case, the model is modified to calculate the state errors between the or- dered and the present value. The new error model with proportional gain is ẋ = Kv(edis cos(✓)� x(t)) ẏ = Kv(edis sin(✓)� y(t)) ˙✓ = Kh(eang � ✓(t)) edis = q (xg � x(t))2 + (yg � y(t))2 eang = arctan ✓ 2 yg � y(t) xg � x(t) ◆ , (2) where Kv is the velocity control gain, Kh is the rotational velocity control gain, edis is the distance error, eang is the 92 angular error, and xg, yg are the Cartesian coordinates to reach. Quasi-Polynomial Representation of Nonlinear Systems Some basic notions of quasi-polynomial and Lotka- Volterra systems are summarised in this section. Generalised Lotka-Volterra Form Representing an ODE in generalised Lotka-Volterra (GLV) form can increase the structural simplicity in ex- change for increasing its dimension. The GLV or quasi- polynomial (QP) form: ẋ i = x i ✓ � i + mP i=1 A ij · nQ k=1 x Bij k ◆ , i = 1, . . . , n, m � n (3) where A and B are n ⇥ m, m ⇥ n real matrices, and � 2 Rn is a vector. The set of non-linear ODEs can be embedded into QP form if it meets two requirements: (i) The non-linear ODEs should follow this form: ẋ s = P is1...isn,js a is1...isn,jsx is1 1 . . . xisn n f(x̄)js , x s (t0) = x0 s , s = 1 . . . n (4) where a is1...isn,js✏R, s = 1 . . . n, and f(x̄) is some scalar function which cannot be reduced to quasy-monomial form. (ii) The partial derivatives of the system Eq.(4) should fulfil: @f @x s = P es1...esn,es b es1...esn,esx es1 1 . . . xesn n f(x̄)es , (5) where b es1...esn,es✏R, s = 1 . . . n. By embedding, we in- troduce the new auxiliary variable: y = fq nQ s=1 xps s . q 6= 0 (6) Differentiating the new, substituted equations we get the QP representation of the original equation Eq.(4): ẋ s = x s 2 4 X is1...isn,js ⇣ a is1...isnjsy js/q nY k=1 xisk��sk�jspk/q !# , s = 1 . . . n (7) where � sk = 1 if s = k and 0 otherwise. A new additional dimension appears as the ODE of the new variable y: ẏ = y " X is1...isn,js ✓ p s x�1 s ẋ s + + X is↵,js es↵,es a is↵,jsbis↵,jsqy (es+js�1)/q ⇥ ⇥ nY k=1 x isk+esk+(1�es�js)pk k ◆# . ↵ = 1 . . . n (8) It is important to mention that the new ODE is not unique because we can choose the parameters p s and q. The quasi-monomial transformation is defined as: x i = nQ k=1 xCik k , i = 1 . . . n (9) where C is an arbitrary invertible matrix. The matrices of GLV can be modified as ˆB = B · C, ˆA = C�1 · A, and ˆ� = C�1 ·�, and the transformed set is also in GLV form. Lotka-Volterra Models The above family of models is split into classes of equiva- lence according to the values of the products M = BA and N = BL. The Lotka-Volterra form gives the rep- resentative elements of these classes of equivalence. If rank(B) = n, then the set of ODEs in Eq.(3) can be embedded into the following m-dimensional set of equa- tions, the so-called Lotka-Volterra model: ż j = z j N j + mX i=1 M ji z i ! , j = 1, . . . ,m (10) where M = BA, N = BL, and each z j represents a so-called quasi-monomial: z j = nY k=1 y Bjk k , j = 1, . . . ,m. (11) Input-Affine QP Model The well known input-affine model of nonlinear sys- tems is given in the following state space model: ẋ = f(x)+ pP j=1 gj(x)uj where f 2 Rn ! Rn, and gj 2 Rn ! Rn are QP functions and the input variable u is p-dimensional. In 93 Figure 1: Membership functions for angular error this case we need a QP form with the matrices from the GLV form ODE set: ẋ i = x i 0 @� i + mX j=1 A0 ij nY k=1 x Bjk k 1 A + + pX l=1 x i 0 @µ li + mX j=1 C lij nY k=1 x Bjk k 1 Au l . i = 1 . . . n (12) It can be proven that if u = H(x) state feedback is in QP form, the closed-loop system remains in QP form as well, but the quasi-monomials of the system will be greater, than the system’s without the feedback. Global Stability Analysis Global equilibrium points can be obtained by finding a Lyapunov function V (·). For LV systems there is a well known Lyapunov function family: V (z) = mX i=1 c i ✓ z i � z⇤ i � ln z i z⇤ i ◆ , c i � 0, i = 1 . . .m (13) and the time derivative: ˙V z = @V (z) @z · @z @t = = 1 2 (z � z ⇤ )(C M +M T C)(z � z ⇤ ), (14) where z ⇤ = (z⇤1 . . . z ⇤ m ) T is the unique positive equilib- rium point, C = diag(c1, . . . , cm) and M is the invari- ant coefficient matrix of the LV form. If CM T +M C is negative semi-definite then z ⇤ is stable and if nega- tive definite then z ⇤ is asymptotically stable. The global stability analysis is thus equivalent to the linear matrix inequality CM T +M C  0 C > 0. (15) The presented Lyapunov function Eq.(13) can be ex- tended to GLV systems by embedding the system using the LV coefficient matrix M = BA. It is necessary to solve the LMI system Eq.(15) for the stability analysis of the QP and LV system: ✓ C 0 0 �CM T �M C ◆ > 0 (16) Controller Design in Quasi-Polynomial Representation The globally stabilizing QP state feedback design prob- lem for QP systems can be formulated as follows [7]. Consider arbitrary quasi-polynomial inputs in the form: u l = rX i=1 k il q̂ i , l = 1 . . . , p (17) where q̂ i = q̂ i (y1, . . . , yn), i = 1, . . . , r are arbitrary quasi-monomial functions of the state variables of the system and k il is the constant gain of the quasi-monomial function q̂ i in the l-th input u l . The closed-loop system will also be a QP system furthermore, the closed-loop LV coefficient matrix M̂ can also be expressed in the form M̂ = B̂  = M0 + pX l=1 rX i=1 k il Mil. (18) Then the global stability analysis of the closed-loop sys- tem with unknown feedback gains k il leads to the follow- ing bilinear matrix inequality M̂ T C +CM̂ = M0 T C +CM0+ pX l=1 rX i=1 k il ⇣ Mil T C +CMil ⌘  0. (19) The variables of the BMI are the p ⇥ r k il feedback gain parameters and the c j , j = 1, ..,m parameters of the Lyapunov function. If the BMI above is feasible then there exists a globally stabilizing feedback with the se- lected structure. Quasi-Polynomial Control of a Differential Drive Mobile Robot Quasi-Polynomial Representation of the Kinematic Model With the given differential drive robot model Eq.(1) and the explained error model Eq.(2) the QP representation can be built. The fist step as mentioned, is to find the new auxiliary variables, that help to eliminate the non-QP ex- pressions. The newly chosen auxiliary variables are: ↵ = cos(✓), � = sin(✓), � = arctan ✓ 2 yg � y xg � x ◆ , � = (xg � x)2 + (yg � y)2, and ✏ = p � + (xg � x). (20) 94 �2 0 2 4 6 0 2 4 6 start finish x [m] y [ m ] PID QP Fuzzy Figure 2: Comparison of different controllers By substituting them into the original system Eq.(2) and performing the differentiation of the new equations Eq.(20): ẋ = x Kv p �↵ x �Kv ! ẏ = y Kv p �� y �Kv ! ˙✓ = ✓ ✓ Kh� ✓ �Kh ◆ ↵̇ = ↵ ✓ �� ↵ ◆ ˙� = � ✓ ↵ � ◆ (21) �̇ = �· ⇣ � 2Kv� p � �✏ + Kvy �✏ �Kvxg� �✏ + Kvx� �✏ + Kvxgy �✏ p � �Kvxy �✏ p � + Kvyg↵ �✏ �Kvy↵ �✏ �Kvygx �✏ p � + Kvxy �✏ p � ⌘ (22) ˙� = � ⇣ � 2Kvxg↵p � + 2Kvxgx � + 2Kvx↵p � � 2Kvx 2 � � 2Kvyg�p � + 2Kvygy � + 2Kvy�p � � 2Kvy 2 � ⌘ (23) ✏̇ = ✏ ⇣ �Kvxg↵ ✏ + Kvx↵ ✏ + Kvxgx ✏ p � �Kvx 2 ✏ p � �Kvyg� ✏ + Kvy� ✏ + Kvygy ✏ p � �Kvy 2 ✏ p � + Kv↵ ✏ �Kvx ✏ p � ⌘ (24) Figure 3: Membership functions for distance error The closed-loop QP system Eq.(21-24) has 29 quasi- monomials. x : ↵ p �x�1 y : � p �y�1 ✓ : �✓�1 ↵ : �↵�1 � : ↵��1 � : � p ���1✏�1, y��1✏�1, ���1✏�1, x���1✏�1, y��1��1/2✏�1, ↵��1✏�1, y↵��1✏�1, x��1��1/2✏�1 � : ↵��1/2, x��1, x↵��1/2, x2��1, ���1/2, y��1, y���1/2, y2��1 ✏ : ↵✏�1, x↵✏�1, x✏�1��1/2, x2✏�1��1/2, �✏�1, y�✏�1, y✏�1��1/2, y2✏�1��1/2, (25) The QP system Eq.(21-24) can obtain the GLV form in- variants with the help of monomials. These matrices are: A the coefficient matrix, B the exponential matrix and � the constant’s matrix. The in- variant matrix product BA results in the LV coefficient matrix M , which is necessary for the global stability analysis. Controller Design As the original model Eq.(2) already contains the propor- tional controller parameters Kh and Kv, they naturally appear in the QP and Lotka-Volterra forms, respectively. The closed-loop Lotka-Volterra coefficient matrix M contains the proportional gains in a linear manner Eq.(26) [7] M = M0 +Kh Mh +Kv Mv, (26) where M 2 R29⇥29. This means, that the globally sta- bilizing feedback design BMI can be formulated, and checked for feasibility. The controller design is practically solved as a control-Lyapunov function-based state feedback. The prescribed Lyapunov function parameters are c1,2,3 = 0.1, c4,5 = 1 and c6�29 = 0. The control gains are ob- tained by solving an LMI version of the feedback design BMI by substituting C into Eq.(19). The obtained values are Kh = 2.78 and Kv = 0.8. 95 Table 1: Rules for the fuzzy controller ^ ^ eang edis !left !right ZERO ZERO ) ZERO ZERO POS ZERO POS NEG ZERO NEG POS ZERO POS POS POS POS POS POS ZERO POS ZERO POS Comparison of Results for Different Controller Types The first and most basic method was the proportional controller approach. This was explained in the second chapter. This is a more systematic approach, but with the increasing non-linearities it is becoming a more difficult method to solve. Fig.3 shows that the robot moves on a significantly wider curve compared to the other two de- signs. Secondly a fuzzy controller was designed to com- pare the original design to. This is a much more intuitive approach, the physical nonlinearities were easily imple- mentable with the cost of inaccuracy. This can be seen in Fig.3, where the robot travels on a much smaller curve. Table 1 shows the definition of the rules. The member- ship functions follow simple triangular shapes according to the acceptable rate of angular and velocity changes. The QP controller performance can be seen in Fig.3. where a classical PID controller and a fuzzy logic con- troller are also shown. It is apparent, that the QP con- troller performs well by means of the length of the trajec- tory. Conclusions A quasi-polynomial representation-based nonlinear con- trol design method has been applied to the kinematics of a differential steering mobile robot in this work. It has been shown, that the quasi-polynomial model of the mo- bile robot’s kinematics extended with the tracking error dynamics has 29 quasi-monomials, i.e. the sizes of matri- ces appearing in the bilinear matrix inequality that should be solved for a globally stabilizing feedback controller Acknowledgement This research was supported by the European Union and the State of Hungary, co-financed by the European Social are reasonably great. The resulting controller and the tra- jectory were compared to the trajectories of a reference PID controller and a fuzzy logic controller. It can be seen that the controller performs satisfactorily. Possible future directions include the application of the LMI and /or BMI cost function to introduce some optimality measure to the problem to be minimised. A next step would be the appli- cation of a time dependent goal position (i.e. a trajectory). Fund in the framework of the TÁMOP 4.2.4. A/2-11-1- 2012-0001 ’National Excellence Programme’. REFERENCES [1] HUANG J., WEN C., WANG W., JIANG Z.-P.: Adap- tive output feedback tracking control of a nonholo- nomic mobile robot, Automatica, 2004, 50(3), 821– 831 [2] MOHARERI O., DHAOUADI R., RAD A.B.: Indirect adaptive tracking control of a nonholonomic mobile robot via neural networks, Neurocomputing, 2012, 88, 54–66 [3] LI J., GUO X., LI Z., CHEN W.: Stochastic Adap- tive Optimal Control of Under-actuated Robots Us- ing Neural Networks, Neurocomputing, 2014, 142, 190-200 [4] HERNÁNDEZ-BERMEJO B., FAIRÉN V.: Nonpoly- nomial vector fields under the Lotka-Volterra normal form, Phys. Lett. A, 1995, 206, 31–37 [5] MAGYAR A., FODOR A.: Quasi-polynomial Control of a Synchronous Generator, Hung. J. Ind. Chem., 2013, 41(1), 55–61 [6] TZAFESTAS S.G.: Mobile Robot Control I: The Lyapunov-Based Method, in Introduction to Mobile Robot Control (TZAFESTAS S.G., ed.) Elsevier, Ox- ford, 2014, 137–183 [7] MAGYAR A., SZEDERKÉNYI G., HANGOS K. M.: Globally stabilizing feedback control of process systems in generalized Lotka-Volterra form, J. Proc. Contr., 2008, 18(1), 80–91 Advertise upcoming meetings, conferences, workshops, new products or services; make public announcements, introduce your research laboratory on the pages of the Hungarian Journal of Industry and Chemistry Please contact us if interested! EDITORIAL OFFICE: UNIVERSITY OF PANNONIA P.O. BOX 158, VESZPRÉM H-8201 (HUNGARY) Tel.: +36 (88) 624-746, E-mail: hjic@almos.uni-pannon.hu; web: hjic.mk.uni-pannon.hu Felelős szerkesztő: Szilágyi Róbert Károly, PhD Kiadja: Pannon Egyetem, 8200 Veszprém, Egyetem u. 10. Levélcím: H-8201 Veszprém, Postafiók 158, Tel.: (88) 624-000 Felelős kiadó: a Pannon Egyetem, Mérnöki Kar dékánja HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 42(2) pp. 97–101 (2014) A NETWORK MODEL FOR SIMULATING THE DYNAMIC BEHAVIOUR OF AN ENERGY DISTRIBUTION NETWORK LEVENTE BÁLINT,1! TIBOR DULAI,2 ÁGNES STARK-WERNER,2 GERGELY APOSTOL,1 AND ÁDÁM TOMPOS1 1 Faculty of Information Technology, University of Pannonia, Egyetem u. 10., Veszprém, 8200 HUNGARY !E-mail: balintlevente.z@gmail.com 2 Department of Electrical Engineering and Information Systems, Faculty of Information Technology, University of Pannonia, Egyetem u. 10., Veszprém, 8200 HUNGARY The aim of this paper is to present the planning and implementation of a network model for an energy distribution network for simulating its dynamical behaviour. The model includes different types of energy sources, distribution centres and consumers. The prepared simulation environment makes it possible to add/remove/modify any type of node – to/form/in the model. Other actions that may be carried out are the initiation/termination of an energy source or change to its production level. The main goal of the simulation is to investigate the effect of an immediate event (such as a power plant failure, decreased production, etc.) on the network and search for the most appropriate substitute possible. The search may happen either on cost bases or time bases. The system created is based on web technology, any user is able to create/save/load his own model and do simulations on it. This work serves as a basis for further work, which involves more complex electrical knowledge and possibilities regarding the network. Keywords: energy distribution network, energy sources, dynamical behaviour, network model, simulation Introduction Energy efficiency has an increasing importance due to cost- and environment-related factors. At the same time, the customer demands have to be satisfied. These circumstances result in a largely heterogeneous energy producing and distributing environment. Renewable energy is applied to a considerable degree in addition to the electricity produced by traditional power plants [1]. This resulted in a system, where different types of energy sources – with different properties, like cost, response time, etc. – are applied in a mixed way. It is a non-trivial task to decide, which sources to apply even in the case of pre-determined demands. Moreover, in practice the amount of the demand appears only when needed. Earlier efforts implemented simulation environments for analysing the effect of the hybrid manner of an electric network [2, 3]. However, the problem becomes more complex when an immediate negative event happens, e.g. an energy source becomes unavailable or has to decrease its output suddenly. The optimal answer depends on many factors: e.g. how quickly a new energy source can be applied or how much it will cost. Moreover, some types of energy sources are only able to produce discrete levels of output power. This paper describes a web-based system, which makes it possible to create and modify an energy distribution network (including different types of energy sources and customers, too) with the goal of simulating different immediate events. The requirements to be collected and the network model and its elements are introduced, followed by details of the implementation and simulation. Requirements: The Desired Properties of the System Functional requirements include mainly the creation and manipulation of network models and simulations. Fig.1 presents the possible applications of the system. The interface contains a tree graph-based model whose root is the power source, the inner nodes are the network stations and the leaves are the consumers. The graph’s edges are the connections between the components. There are multiple projects in the system, each project contains a model. The main functions of the system are available through the modeling and simulation interfaces. During the simulation the components can be switched on and off. After two hours of simulation a statistic diagram can be viewed for each component of the model, showing its consumption or production. Information about components can be edited while the system is not simulating. Next to the functional requirements, the system-related requirements are that: • the system has to be web-based, • no user authentication is needed, • the reaction time is less than 1 s, even if there are 100 users, and • any system or communication failure has to be presented by an unambiguous alert to the user. 98 The Interface The interface of the system contains four important elements. On the upper edge is the modelling control bar with the component adding buttons, background theme switcher and save/load buttons. By pressing the load button, the system opens the project selector. The save button must be pressed before loading another project in order to save the current one. Under the right side of the control bar, there is the simulation control panel. There are buttons on it for editing corporate and residential priority orders, and sliders for the importance of certain properties while calculating the priority order of the power sources. There are buttons for starting, stopping, and restarting simulation time. The simulation time begins at 6:00 am and the step time can be set from 1 to 60 minutes. In the bottom left-hand corner, the event-log is located, which indicates the basic events like start/stop simulation, add/remove components in black and important events like system failure or insufficient power in red. All events have a time-stamp. While hovering the mouse over a component, the system shows the detailed properties of the element. Double clicking selects a component or deletes a connection. While a component is selected, the edit, delete and diagram buttons appear in the modelling control bar. The diagram button works only after a minimum of two hours of simulation. In simulation mode the modelling control bar is hidden and double clicking switches a component on/off. The Network Model and its Elements As the system’s main potential is in handling different types of network nodes and making the investigation of a heterogeneous energy producing and distribution system possible, the elements of the network models and their properties have a huge role. The network nodes can be classified as follows. Energy Sources The system considers renewable energy sources (water, bio, wind, solar, and geo), non-renewable energy sources (natural gas, petrol, coal, nuclear), and quick start-up energy sources (gas turbines) when the energy provided by power plants is insufficient for consumers. Each energy source has a unique name and a maximum output, which defines an upper limit for the source (Fig.2). A minimum output is defined as a lower limit, below that the operation is not profitable. Start-up cost is the cost needed to switch on the source, and start- time is the time required to start producing energy. The operational cost is a financial value which is needed for one hour of operation. Power sources have an environmental multiplier. This multiplier is 1 for renewable sources, 0.2 for non-renewable ones and 0.1 for quick start-up sources by default, but it can be edited any time during the modelling. Larger values are better, the maximum is 1 but the value cannot be 0. The system calculates an order between sources from their start-up cost, operational cost and environmental multiplier. This order can be modified by the user before starting the simulation. Distributing Elements Distributing elements for the system are power plant stations that distribute the energy to multiple network stations that further distribute the energy to retail stations directly transmit the energy to consumers. Distributing stations have only two properties (Fig.3). They have a unique name to identify them and a maximum throughput, which is the maximum power that the station can distribute to lower stations or consumers. If that limit is reached the station will not pick up more power even if the consumers need more. Figure 1: Example uses of the prepared system Figure 3: Edit station properties (Név: name or location, Teljesítmény: power) Figure 2: Edit power source properties (Név: name or location, Bekapcsolási költség: cost of starting up, Indítási idő: start-up time, Üzemköltség: operational cost, Környezeti szorzó: environmental factor, Optimális sorrend: optimal order) 99 Consumers The residential consumers are linked to the model through homes, while corporate consumers include factories, public institutions, etc. Consumers have a rate of consumption that describes their maximum required power. Current consumption is estimated during the daytime in the simulation. Consumers also have an order of priority that can be modified by the user but corporate consumers always have priority over residential consumers (Fig.4). Connections Connections have a length property, which will be used for calculating the loss, dissipated by the cables resistance. Networks can be built up by the previously introduced network nodes by connecting them appropriately. The order of connectivity is as follows: • power sources, • power plant distribution station, • network distribution station, • retail distribution station, and • consumer. Connections can be established only in this order regardless of the performance of individual components. Database and Implementation Details Properties of the network model elements and the network formed are stored in an associative array that can be used directly by the system. At each time it is saved, this array is stored in JavaScript Object Notation (JSON) format in a text file for each project separately. The database array contains an entry for all components of the model. Every component has an ID, a name, a type and position on the modelled field. Energy sources also have a fuel type, start-up time, start-up cost, running cost, environmental factor, and minimum and maximum performance level attributes. Stations have a defined maximum tolerance. For consumers, there is a stored average consumption rate, which dynamically changes over time. Each object has a sub-array, which stores the ID of any other components that are connected to the current one. Consumer objects have priority numbers that define the relevance of the consumer. If the power source cannot produce enough energy to satisfy all needs, then the consumer that has the highest priority number will be switched off from the network. The priority order equals the order of creation of components, but the corporate consumers always have preference. The order can be changed by the user any time in modelling mode and during the simulation. Implementation Details The program is implemented in JavaScript language with plug-ins like jQuery and jsPlumb. The plug-in jQuery is the most popular extension of the JavaScript language. It makes it easier to select, move, remove or modify the elements on the interface by selecting the element by its ID and modifying its properties. The system uses an in-house jQuery plug-in, windows.js. This plug-in creates free movable small windows on the interface and fills them with appropriate content defined by the programmer. The jsPlumb plug-in is an open- source plug-in for JavaScript, created to visually connect elements on the interface with straight lines or Bezier curves. The database is stored in JSON format text files, created and edited by Hypertext Preprocessor (PHP) functions and called asynchronously by JavaScript’s Asynchronous JavaScript And XML (AJAX) protocol. PHP is a server-side scripting language, but in our case its only function is to properly save and load database files to the server computer. The whole software can be started in a browser on the user’s computer and it connects to the server only while loading or saving the project. The interface is fully interactive, in modelling mode the objects can be moved or removed, in simulating mode switched on and off. Simulation Possibilities The system was created for simulating immediate events during the operation of a modelled energy distribution network. After forming the network, the software makes it possible to cause different immediate events. These severe a connection or switch on/off sources, stations and/or consumers. All of these events can be carried out on the graphical interface of the software (Fig.5). The simulation is daytime-based. It starts at 6:00 am. Users can set the step time from 1 min or sec to 60 min or sec. At the beginning the system distributes the energy produced by the power plant(s) to all consumers in a user defined order. If a power plant does not produce enough power then the last consumer will be switched off and so on. The electricity-related behaviour will be implemented in the future into the software. A possible solution is to build an ideal flow model as a cooperative game over a graph with the sources and consumers located at the nodes, each described by a maximum supply or desired demand and the power lines Figure 4: Edit consumer properties (Név: name or location, Fogyasztás: consumption, Prioritás: priority) 100 represented by the edges, each with a given power transmission capacity and admittance value describing their ability to transmit electricity. [4] The user is able to initiate immediate events and the system reacts to them in real-time and recalculates the distribution. Calculation of the consumption rate is based on the Hungarian average with a small random factor. An Example for a Simulation Let us assume that we have a power plant with a maximum power of 100 kW. The model has one of each type of stations, three households (10-10 kW maximum) and a factory (80 kW maximum) (Fig.6). Normally, all of the consumers are satisfied, but in the peak times, when the factory and the houses also approach their maximum consumption, the power will not be enough for every consumer, so the system switches off the retail consumer that is last in the priority order. There will not be a consumption reduction, but a complete shutdown. Subsequently, the plant will have to produce 95% of the maximum power. In this example all of the stations can transfer the maximum of 100 kW, but if a station cannot transfer this, then the consumer with the largest priority number connected to this station will be switched off, regardless of the current level of the power source. If the consumption is too high, it is possible that complete distribution stations are switched off. Further Developments The current system can model and simulate a fictional power network with components and connections, but the main goal of the project is to create a simulating software that can model and simulate a real electrical network. It is important to emphasize that the unit of the power of the sources and the consumption are now only kW. We need to include that the voltage of connection cables depends on the length of the connection and the diameter of the cable that are required to calculate the amperage from the resistance of the cable. International connections transmit 750 kV, and the main network in Hungary transmits 400 kV starting from the power plants (Fig.7). The distribution stations transform and divide it to lower values until the voltage reaches the necessary 250 V for retail and 360 V (or more) for corporate consumers. During the transformations, the system needs to use Kirchhoff’s current law that states that at any node in an electrical circuit, the sum of currents flowing into that node is equal to the sum of currents flowing out of that node. That means that the system has to calculate amperages for every incoming and outgoing connection of the distribution stations based on the consumption of the consumers and the target voltage of the cables. The biggest problem of power distribution networks in reality is to satisfy all consumers at all the time, even if the load is very low, without making large changes to the output of the power sources everyday. The current solution in Hungary is nighttime power. Households using electric water heaters have nighttime electric meters, when the consumption of corporate consumers is drastically reduced. This solution slightly optimizes the operation of power plants. Figure 5: The graphical interface of the software (Új erőmű: new powerplant, Új fogyasztó: new consumer, Új állomás: new station, Új kapcsolat: new connection, Mentés: Save, Betöltés: Load, Lakossági sorrend: order of households, Céges sorrend: order of companies, MaxKöltség: maximum cost, Környezet: environment, Bekaps: start-up, Indulás: begin, Indítás: start, Leállítás: stop, Újraindítás: restart Figure 6: The model of our example (méter: distance in m) 101 Conclusion This paper presented the planning phase and implementation of a web-based system that makes it possible to construct an energy distribution network together with different types of energy sources and consumers. The simulations are carried out to determine the most effective response to different immediate events (e.g. a sudden network problem or failure of an energy source). The simulation takes into account both cost-related and time-related properties, too. The present state of the software provides a basis for further works, which intend to add deeper electricity-related properties and behaviour to the system. ACKNOWLEDGEMENT The authors are thankful for scientific guidance by Prof. KATALIN M. HANGOS from the Computer and Automation Research Institute. We acknowledge the financial support of the Hungarian State and the European Union under the TAMOP-4.2.2.A- 11/1/KONV-2012-0072 project. REFERENCES [1] BREITNER M.H.: Modeling the transformation of the German energy system until 2050: A multi-criteria, long term optimization problem with many constraints, Operations Res. Proc. 2012, 369-374 [2] GEYSEN D., BOOIY P., WARMER C.: A framework for simulation and control of hybrid energy networks, Proc. EnergyCon 2014, Dubrovnik, Croatia, 2014 [3] PADULLES J., AULT G.W., MCDONALD J.R.: An integrated SOFC plant dynamic model for power systems simulation, Int. J. Power Sources 2000, 86, 495-500 [4] CSERCSIK D., KÓCZY L.Á.: Externalities in the games over electrical power transmission networks, IEHAS Discussion Papers, 2011, No 1125 Figure 7: The Hungarian transmission network Advertise upcoming meetings, conferences, workshops, new products or services; make public announcements, introduce your research laboratory on the pages of the Hungarian Journal of Industry and Chemistry Please contact us if interested! EDITORIAL OFFICE: UNIVERSITY OF PANNONIA P.O. BOX 158, VESZPRÉM H-8201 (HUNGARY) Tel.: +36 (88) 624-746, E-mail: hjic@almos.uni-pannon.hu; web: hjic.mk.uni-pannon.hu Felelős szerkesztő: Szilágyi Róbert Károly, PhD Kiadja: Pannon Egyetem, 8200 Veszprém, Egyetem u. 10. Levélcím: H-8201 Veszprém, Postafiók 158, Tel.: (88) 624-000 Felelős kiadó: a Pannon Egyetem, Mérnöki Kar dékánja HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 42(2) pp. 103–107 (2014) ON THE PARAMETRIC UNCERTAINTY OF WEAKLY REVERSIBLE REALIZATIONS OF KINETIC SYSTEMS GYÖRGY LIPTÁK B1, GÁBOR SZEDERKÉNYI1,2 AND KATALIN M. HANGOS1,3 1Process Control Research Group, MTA SZTAKI, Kende u. 13-17, Budapest, 1111, HUNGARY 2Faculty of Information Technology, Péter Pázmány Catholic University, Práter u. 50/a, Budapest, 1083, HUNGARY 3Department of Electrical Engineering and Information Systems, University of Pannonia, Egyetem u. 10, Veszprém, 8200, HUNGARY BE-mail: lipgyorgy@gmail.com The existence of weakly reversible realizations within a given convex domain is investigated. It is shown that the domain of weakly reversible realizations is convex in the parameter space. A LP-based method of testing if every element of a convex domain admits weakly reversible realizations is proposed. A linear programming method is also presented to compute a stabilizing kinetic feedback controller for polynomial systems with parametric uncertainty. The proposed methods are illustrated using simple examples. Keywords: parametric uncertainty; computational methods; optimization; kinetic systems Introduction The notion of parametric robustness is well-known and central in linear and nonlinear systems and control theory [1]. It is used for ensuring a desirable property, such as stability, in a given domain in the parameter space around a nominal realization having the desired property. The aim of the paper is to extend the notions and tools of parametric robustness for a class of positive polyno- mial systems, namely a class of kinetic systems. Only the very first steps are reported here that offer a computation- ally efficient method for checking one of the many impor- tant properties of kinetic systems, their weak reversibility. Basic Notions and Methods The basic notions and tools related to reaction kinetic sys- tems and their realizations are briefly summarized in this section. Kinetic Systems, their Dynamics and Structure Deterministic kinetic systems with mass action kinetics or simply chemical reaction networks (CRNs) form a wide class of non-negative polynomial systems, that are able to produce all the important qualitative phenomena (e.g. stable/unstable equilibria, oscillations, limit cycles, multiplicity of equilibrium points and even chaotic be- haviour) present in the dynamics of nonlinear processes [2]. The general form of dynamic models studied in this paper is the following ẋ = M · (x), (1) where x 2 Rn is the state variable and M 2 Rn⇥m. The monomial vector function : Rn ! Rm is defined as j(x) = nY i=1 x Yij i , j = 1, . . . ,m (2) where Y 2 Nn⇥m 0 . The system Eq.(1) is kinetic if and only if the matrix M has a factorization M = Y ·Ak. (3) The Kirchhoff-matrix Ak has non-positive diagonal and non-negative off-diagonal elements and zero column sums. The matrix pair (Y,Ak) is called the realization of the system Eq.(1). The chemically originated notions: The chemically originated notions of kinetic systems are as follows: the species of the system are denoted by X1, . . . , Xn, and the concentrations of the species are the state variables of Eq.(1), i.e. xi = [Xi] � 0 for i = 1, . . . , n. The structure of kinetic systems is given in terms of its com- plexes Ci, i = 1, . . . ,m that are non-negative linear combinations of the species i.e. Ci = Pn j=1[Y ]jiXj for i = 1, . . . ,m, and therefore Y is also called the complex composition matrix. The reaction graph: The weighted directed graph (or reaction graph) of kinetic systems is G = (V,E), where 104 V = {C1, C2, . . . , Cm} and E denote the set of ver- tices and directed edges, respectively. The directed edge (Ci, Cj) (also denoted by Ci ! Cj) belongs to the re- action graph if and only if [Ak]j,i > 0. In this case, the weight assigned to the directed edge is Ci ! Cj is [Ak]j,i. Stoichiometric subspace: Stoichiometric subspace S is given by the span of the reaction vectors S = {[Y ]·i � [Y ]·j | [Ak]ij > 0}. (4) The stoichiometric compatibility classes of a kinetic sys- tem are the affine translations of the stoichiometric sub- space: (x0 + S) \ Rn �0. Structural Properties and Dynamical Behaviour It is possible to utilize certain structural properties of ki- netic systems that enable us to effectively analyze the sta- bility of the system. Deficiency: There are several equivalent ways to de- fine deficiency. We will use the following definition � = dim(ker(Y ) \ Im(BG)), (5) where BG is the incidence matrix of the reaction graph. It is easy to see that deficiency is zero if ker(Y ) = {0} or equivalently rank(Y ) = m. Weak reversibility: A CRN is called weakly reversible if whenever there exists a directed path from Ci to Cj in its reaction graph, then there exists a directed path from Cj to Ci. In graph theoretic terms, this means that all components of the reaction graph are strongly connected components. Deficiency zero theorem: A weakly reversible kinetic system with zero deficiency has precisely one equilibrium point in each positive stoichiometric compatibility class that is locally asymptotically stable (conjecture: globally asymptotically stable). Computing Weakly Reversible Realizations Formulated as an Optimization Problem In this section, first a method for computing weakly re- versible realization based on Ref.[3] is briefly presented. We assume that we have a kinetic polynomial system of the form Eq.(1). We use the fact known from the literature that a real- ization of a CRN is weakly reversible if and only if there exists a vector with strictly positive elements in the kernel of Ak, i.e. there exists b 2 Rn + such that Ak · b = 0 [4]. Since b is unknown, too, this condition in this form is not linear. Therefore, we introduce a scaled matrix Ãk Ãk = Ak · diag(b) (6) where diag(b) is a diagonal matrix with elements of b. It is clear from Eq.(6) that Ãk is also a Kirchhoff matrix and that 1 2 Rm (the m-dimensional vector containing only ones) lies in kernel of Ãk. Moreover, it is easy to see that Ãk defines a weakly reversible network if and only if Ak corresponds to a weakly reversible network. Then, the weak reversibility and the Kirchhoff property of Ãk can be expressed using the following linear constraints mX i=1 h Ãk i ij = 0, j = 1, . . . ,m mX i=1 h Ãk i ji = 0, j = 1, . . . ,m h Ãk i ij � 0, i, j = 1, . . . ,m, i 6= j h Ãk i ii  0, i = 1, . . . ,m. (7) Moreover, the equation Eq.(3) is transformed by diag(b) (we can do this, because diag(b) is invertible): M · diag(b) = Y ·Ak · diag(b)| {z } Ãk (8) Finally, by choosing an arbitrary linear objective function of the decision variables Ãk and b, weakly reversible real- izations of the studied kinetic system can be computed (if any exist) in a LP framework using the linear constraints Eq.(7) and (8). Weakly Reversible CRN Realizations In this section, first the convexity of the weakly reversible Kirchhoff matrix will be shown. After that the practical benefits of this property will be demonstrated in the field of system analysis and robust feedback design. Convexity of the Weak Reversibility in the Parameter Space Theorem 1. Let A (1) k and A (2) k be m ⇥ m weakly re- versible Kirchhoff matrices. Then the convex combination of the two matrices remains weakly reversible. Proof. The idea behind the proof is based on Ref.[5]. A Kirchhoff matrix is weakly reversible if and only if there is a strictly positive vector in its kernel. Therefore strictly positive vectors p1, p2 exist such as A(1) k · p1 = 0 and A (2) k · p2 = 0. Let us define the following scaled Kirchhoff matrix: Â(1) k = A (1) k · diag(p1) and  (2) k = A (2) k · diag(p2). These scaled matrices have identical structures to the original ones. Moreover, Â(1) k ·1(m) = 0 and  (2) k ·1(m) = 0 where the vector 1(m) denotes the m dimensional column vector composed of ones. For that (�Â(1) k + (1� �)Â(2) k ) · 1(m) = 0. (9) for any � 2 [0, 1]. Therefore the convex combination of the original two realizations has to be weakly re- versible. 105 (a) (b) (c) Figure 1: A weakly reversible reaction graphs of the three realizations (Y,A(1) k ), (Y,A(2) k ) and (Y,A(3) k ) Weak Reversibility of CRN Realizations with Parametric Uncertainty We assume that a CRN with parametric uncertainty is given as ẋ = M · (x), (10) where x 2 Rn is the state variable, 2 Rn ! Rm contains the monomials and the matrix M 2 Rn⇥m is an element of the following set M = � lX i=1 ↵iMi | (8i : ↵i � 0) ^ lX i=1 ↵i = 1 . (11) The goal is to find a method for checking the weak re- versibility of the system Eq.(10) for all matrices M 2 M. When all vertices Mi have a weakly reversible real- ization (Y,A(i) k ) then any element of the set M has a re- alization (Y,Ak) such that Ak is the convex combination of the Kirchhoff matrices A (i) k . The obtained realization Ak will be weakly reversible due to Theorem 1. There- fore, it is enough to compute a weakly reversible realiza- tion for each matrix Mi by using the previously presented LP-based method. A Simple Example Let us consider the following polynomial system  ẋ1 ẋ2 � = M · 2 4 x1 x2 x1x2 3 5 , (12) where M is an arbitrary convex combination of the fol- lowing three matrices M1 =  0 1 �1 1 �1 0 � , M2 =  �1 1 0 1 �1 0 � , and M3 =  0 1 �1 0 0 0 � . Figure 2: A weakly reversible realization of the convex combination M = 0.2M1 + 0.4M2 + 0.4M3 In order to show weak reversibility for all possible con- vex combinations, we have to find a weakly reversible re- alization for each matrix M1, M2 and M3. The resulting weakly reversible reaction graphs are depicted in Fig.1, while Fig.2 illustrates an inner point realization which is weakly reversible too. Computing Kinetic Feedback for a Polynomial System with Parametric Uncertainty Besides the possible application of the above described LP-based method for robust stability analysis, it can also be used for stabilizing feedback controller design. For this purpose, a generalized version of our preliminary work on kinetic feedback computation for polynomial systems to achieve weak reversibility and minimal defi- ciency [6] is used here. The Feedback Design Problem We assume that the equation of the open-loop polynomial system with linear constant parameter input structure is given as ẋ = M · (x) +Bu, (13) where x 2 Rn is the state vector, u 2 Rp is the input and 2 Rn ! Rm contains the monomials of the open-loop system. The input matrix is B 2 Rn⇥p, the correspond- ing complex composition matrix is Y with rank m, and M 2 Rn⇥m is an element of the following set M = ( lX i=1 ↵iMi | (8i : ↵i � 0) ^ lX i=1 ↵i = 1 ) . (14) Moreover, a positive vector x 2 Rn >0 being the desired equilibrium point is given as a design parameter. Note that the above polynomial system is not necessarily ki- netic, i.e. not necessarily positive, and may not have a positive equilibrium point at all. 106 The aim of the feedback is to set a region in the state space R ✓ Rn �0 where x is (at least) a locally asymptot- ically stable equilibrium point of the closed-loop system for all M 2 M. For this purpose we are looking for a feedback in the form u = K (x) (15) which transforms the open-loop system into a weakly re- versible kinetic system with zero deficiency for all M 2 M with the given equilibrium point x. Feedback Computation Similarly to the realization computation, the matrix K will be determined by solving an LP problem. The con- vexity result shows that it is enough to compute one weakly reversible realization (Y,A(r) k ) in each vertex Mr to ensure weak reversibility for all possible closed-loop systems. All realizations will have zero deficiency, be- cause of the rank condition rank(Y ) = m [7]. First we note, that the realization (Y,A(r) k ) that corre- sponds to the closed-loop system is Mr +B ·K = Y ·A(r) k . (16) where the matrix A (r) k should be Kirchhoff mX i=1 [Ãk]ij = 0, j = 1, . . . ,m [Ãk]ij � 0, i, j = 1, . . . ,m, i 6= j [Ãk]ii  0, i = 1, . . . ,m. (17) In order to obtain a weakly reversible closed-loop system with an equilibrium point x, the matrix A (r) k should be weakly reversible and has to have the vector (x) in its right kernel, i.e. A (r) k · (x) = 0. (18) Finally, by choosing an arbitrary linear objective func- tion of the decision variables A (1) k , . . . , A (l) k and K, the feedback gain K can be computed (if it exists) in a LP framework using the linear constraints Eqs.(16-18). With the resulting feedback gain K, the point x will be an equilibrium point of all possible closed-loop sys- tems, and x will be locally asymptotically stable in the region S = (x+S)\Rn �0, where S is the stoichiometric subspace of the closed-loop system. Example Let the open-loop system be given as ẋ = M 2 4 x1x2 x2x3 x1 3 5+ 2 4 0 1 0 3 5 u (19) Figure 3: Weakly reversible realization of the closed-loop system, where M = 0.6M1 + 0.2M2 + 0.2M3 where M is an arbitrary convex combination of the fol- lowing three matrices: M1 = 2 4 �1 1 0 2 1 2 1 �1 0 3 5 , M2 = 2 4 0 0 0 1 1 3 0 0 0 3 5 , and M3 = 2 4 0 1 �1 2 0 3 0 �1 1 3 5 . The desired equilibrium point x = [1 1 1]T . We are looking for a feedback law with gain K which transforms the matrices Mi into weakly reversible kinetic systems with the given equilibrium point. By solving the feedback design LP optimization prob- lem using the linear constraints Eqs.(16-18), the com- puted feedback is in the following form: u = ⇥ 2 1 2 ⇤ (x). (20) Fig.3 depicts a weakly reversible realization of the closed-loop system. The obtained closed-loop system in an inner point of the convex set M has the following sto- ichiometric subspace: S = span 0 @ 2 4 1 1 0 3 5� 2 4 0 1 1 3 5 , 2 4 1 1 0 3 5� 2 4 1 0 0 3 5 1 A . (21) Therefore, the equilibrium point x will be asymptotically stable with the region S = (x+S)\Rn �0. Note, that one should choose the initial value of the state variables from S . Fig.4 shows the time dependent behaviour of the closed-loop solutions started from different initial points in S . 107 Figure 4: Time-domain simulation of the closed-loop system Conclusion It is shown in this paper that the domain of weakly re- versible realizations is convex in the parameter space. This property is utilized for developing methods in sys- tem analysis and robust control design. An LP-based op- timization method is proposed for testing if every element of a convex domain given by its extremal matrices admits a weakly reversible realization. An LP-based feedback design method is also proposed that guarantees stability with a desired equilibrium point. The proposed methods are illustrated with simple examples. REFERENCES [1] SILJAK D.: Parameter space methods for robust con- trol design: a guided tour, Aut. Contr., IEEE Trans- actions on, 1989, 34(7), 674–688 [2] ANGELI D.: A tutorial on chemical network dynam- ics, Eur. J. Contr., 2009, 15, 398–406 [3] SZEDERKÉNYI G., HANGOS K. M., TUZA Z.: Find- ing weakly reversible realizations of chemical re- action networks using optimization, MATCH Com- mun. Math. Comp. Chem., 2012, 67, 193–212 [4] FEINBERG M., HORN F.: Chemical mechanism structure and the coincidence of the stoichiometric and kinetic subspaces, Arch. Rational Mechanics and Analysis, 1977, 66(1), 83–97 [5] RUDAN J., SZEDERKÉNYI G., HANGOS K.M., PÉNI T.: Polynomial time algorithms to determine weakly reversible realizations of chemical reaction networks, J. Math. Chem., 2014, 52(5), 1386–1404 [6] LIPTÁK G., SZEDERKÉNYI G., HANGOS K.M.: Ki- netic feedback computation for polynomial systems to achieve weak reversibility and minimal deficiency, 13th Eur. Control Conf. ECC, Strasbourg, France, 2014, 2691–2696 [7] ANGELI D., LENHEE, P.D., SONTAG E.D.: On the structural monotonicity of chemical reaction networks, Proc. 45th IEEE Conf. Decision and Control, 2006, 7–12 Advertise upcoming meetings, conferences, workshops, new products or services; make public announcements, introduce your research laboratory on the pages of the Hungarian Journal of Industry and Chemistry Please contact us if interested! EDITORIAL OFFICE: UNIVERSITY OF PANNONIA P.O. BOX 158, VESZPRÉM H-8201 (HUNGARY) Tel.: +36 (88) 624-746, E-mail: hjic@almos.uni-pannon.hu; web: hjic.mk.uni-pannon.hu Felelős szerkesztő: Szilágyi Róbert Károly, PhD Kiadja: Pannon Egyetem, 8200 Veszprém, Egyetem u. 10. Levélcím: H-8201 Veszprém, Postafiók 158, Tel.: (88) 624-000 Felelős kiadó: a Pannon Egyetem, Mérnöki Kar dékánja HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 42(2) pp. 109–113 (2014) STABILITY AND PARAMETER SENSITIVITY ANALYSES OF AN INDUCTION MOTOR ATTILA FODOR1, ROLAND BÁLINTB1, ATTILA MAGYAR1, AND GÁBOR SZEDERKÉNYI2 1Department of Electrical Engineering and Information Systems, University of Pannonia, Egyetem u. 10., Veszprém, 8200, HUNGARY 2Pazmány Péter Catholic University, Práter u. 50/a, Budapest, 1083, HUNGARY BE-mail: balint.roland27@gmail.com A simple dynamical model of an induction motor is derived and analyzed in this paper based on engineering principles that describe the mechanical phenomena together with the electrical model. The used state space model consists of nonlinear state equations. The model has been verified under the usual controlled operating conditions when the speed is controlled. The effect of load on the controlled induction motor has been analyzed by simulation. The sensitivity analysis of the induction motor has been applied to determine the model parameters to be estimated. Keywords: induction motor, stability analysis, sensitivity analysis Introduction Induction motors (IM) are the most commonly used elec- trical rotating machines in several industrial applications. Irrespective of size and the application area, these motors share the most important dynamical properties, and their dynamical models have a similar structure. Because of the specialties and great practical impor- tance of IMs in industrial applications, their modelling for control purposes is well investigated in the literature. Besides basic textbooks [1-3], there are several papers that describe the modelling and use the developed mod- els for the design of different types of controllers: vec- tor control [1, 4], sensorless vector control [5] and direct torque control (DTC) [6]. The aim of this paper is to build a simple dynamical model of the IM and to perform its parameter sensitivity analysis. The results of this analy- sis will be the basis of the next step since the final aim of our study is to estimate the parameters of the IM and design a controller that can control the speed and torque of the IM. The state space model has been implemented in the Matlab/Simulink environment which enables us to analyze the parametric sensitivity based on simulation ex- periments. Nonlinear Model of an Induction Motor In this section a state space model of an induction mo- tor (IM) is presented. The model development is largely based on Refs.[1, 8-10]. For constructing the IM model, the following modelling assumptions are made: 1. a symmetrical triphase stator winding system is as- sumed, 2. the flux density is radial in the air gap, 3. the copper loss and the slots in the machine can be neglected, 4. the spatial distribution of the stator fluxes and aper- tures wave are considered to be sinusoidal, 5. stator and rotor permeability are assumed to be infi- nite with linear magnetic properties. According to the above modelling conditions the mathematical description of the IM is developed through the space-vector theory. If the voltage of the stator is pre- sumed to be the input excitation of the machine, then the spatial distribution along the stator of the x phase volt- age can be described by the complex vector v sx (t). We can determine the orientation of the voltage vector v s , the direction of the respective phase axis and the voltage po- larity. i s (t) = 2 3 � a0 · i sa (t) + a1 · i sb (t) + a2 · i sc (t) � = p 2 · i e↵ (t) · ej!0t+⇡ 2 +' i , (1) where a is the ej120 o vector and i sa , i sb and i sc are the following: i sa (t) = Re(a0 · i s (t)) = Re(i s (t)), i sb (t) = Re(a2 · i s (t)), and i sc (t) = Re(a1 · i s (t)). In Eq.(1), 2/3 is the normalizing factor. The flux den- sity distribution can be obtained by integrating the cur- rent density wave along the cylinder of the stator. The flux linkage wave is a system variable, because it con- tains detailed information about the winding geometry. 110 Figure 1: The equivalent circuit of the IM The rotating flux density wave induces voltages in the in- dividual stator windings. Thus stator voltage v s (t) can be represented as the overall distributed voltages in all phase windings: v s (t) = 2 3 � a0 · v sa (t) + a1 · v sb (t) + a2 · v sc (t) � = p 2 · v e↵ (t) · ej!0t+⇡ 2 +' u , (2) where a is the ej120 o vector and i sa , i sb and i sc are the following: v sa (t) = Re(a0 · u s (t)) = Re(u s (t)), v sb (t) = Re(a2 · u s (t)), and v sc (t) = Re(a1 · u s (t)). Considering the stator of the IM as the primer side of the transformer, then using Kirchoff’s voltage law the follow- ing equation can be written (Fig.1): v s (t) = i s (t) ·R s + d� s (t) dt . (3) As for the secondary side of the transformer, it can be deduced that the same relationship is true for the rotor side space vectors: v r (t) = ir(t) ·Rr + d� r (t) dt = 0. (4) Eqs.(3) and (4) describe the electromagnetic interaction as the connection of first order dynamical subsystems: � s (t) = i s (t) · L s + i r (t) · L m , and (5) � r (t) = i s (t) · L m + i r (t) · L r . (6) Since four complex variables (i s (t) , i r (t) , � s (t) , and � r (t)) are presented in Eqs.(5) and (6), flux equations are needed to complete the relationship between them. The mechanical power (P mech (t)) of the IM can be defined as: P mech (t) = W mech (t) dt , (7) where the mechanical energy (P mech (t)) in the rotating system can be given by the following expression: P mech (t) = W mech (t) dt = T mech (t) · ! r , (8) where T e (t) is the torque and ! r is the angular velocity of the IM. Afterwards the energy balance of the IM is as follows: W e = W mech +W R +W Field , (9) Figure 2: The equivalent circuit of the d-axis of the IM Figure 3: The equivalent circuit of the q-axis of the IM where P e = W e (t) dt = 3 2 Re (u s · i s + u r · i r ) (10) is the input electrical power, P R = W R (t) dt = 3 2 Re � R s · |i s |2 +R r · |i r |2 � (11) represents the resistive power losses and P Field = W Field (t) dt = 3 2 Re ✓ d� s dt i s + d� r dt i r ◆ (12) is the air gap power. Using Eqs.(8-12), P mech (t) = T mech (t) · ! r = 3 2 · Lm L r · � s (t)⇥ i s (t), and (13) T e = 1.5p · (� ds i qs � � qs i ds ). (14) The equivalent circuit of the IM can be decomposed to di- rect axis and quadratic axis components by Park’s trans- formation as shown in Figs.2 and 3. The actual terminal voltage v of the windings can be written in the form v = ± JX j=1 (Rj ij)± JX j=1 ✓ d�j dt ◆ , (15) where ij are the currents, Rj are the winding resistances, and �j are the flux linkages. The positive directions of the stator currents point out of the IM terminals. By composing the d and q axes of the IM the follow- ing equations can be written: v qs = R s · i qs + � qs dt + ! · � ds , (16) v ds = R s · i ds + � ds dt � ! · � qs , (17) v qr = R0 r · i0 ds + �0 dr dt + (! � ! r ) · �0 dr , and (18) v dr = R0 r · i0 dr + �0 dr dt � (! � ! r ) · �0 qr . (19) 111 Figure 4: Response to the step change of the speed-controlled IM d! m dt = 1 2H (T e � F · ! m � T mech ) (20) where ! is the reference frame angular velocity, ! r is the electrical angular velocity, � qs = L s · i qs + L m i0 qr , (21) � ds = L s · i ds + L m i0 dr , (22) � qr = L r · i0 qr + L m i0 qs , and (23) � dr = L r · i0 dr + L m i0 ds . (24) The above model can be written in state space form by ex- pressing the time derivative of the fluxes and ! from the voltage and swing equations Eqs.(21-24). The nonlinear state space model of the IM is given by Eqs.(25-29): d� qs dt = �R s L s � L2 m L0 r · � qs R s · L m L0 r (L s � L2 m L0 r ) · �0 qr � ! · � ds + v qs , (25) d� ds dt = �R s L s � L2 m L0 r · � ds + R s · L m L0 r (L s � L2 m L0 r ) · �0 dr � ! · � qs + v ds , (26) d�0 qr dt = �R0 r L m � L0 r·Ls Lm · � qs + �R0 r · L s L m · (L m � L0 r·Ls Lm ) · �0 qr � ! · �0 dr + ! r · �0 dr + v qr , (27) d�0 dr dt = �R0 r L m � L0 r·Ls Lm · � ds + �R0 r · L s L m · (L m � L0 r·Ls Lm ) · �0 dr + ! · �0 qr � ! r · �0 qr + v dr , and (28) d! r dt = 1 2H · 1.5 · p · L m L m · (L s � L2 m L0 r ) · �0 dr · � qs � 1 2H · 1.5 · p · L m L m · (L s � L2 m L0 r ) · �0 qr · � ds � F 2H · ! r � T m 2H . (29) The state vector of the above model is x=[� qs , � ds , �0 qr , �0 dr , ! r ]T 2 R5, and the input variables are organized in terms of the the input vector u=[v q , v d , �T mech ]T 2 R3. It is assumed that all the state variables can be measured i.e. y = x. Model Validation The dynamical properties of the IM have been investi- gated. The response of the speed-controlled motor has been tested under step-like changes. The simulation re- sults are shown in Fig.4, where the fluxes (� qs , � ds , �0 qr , and �0 dr ) and the angular velocity (!) are shown. Model Analysis The above model Eqs.(21-24) has been verified by simu- lation against engineering intuition. Local Stability Analysis As the final aim of our research is to estimate the param- eters of a particular Grundfos IM, first of all the resis- tances (R s and R r ) of the IM were measured. Afterwards the values of the inductances (L ls , L lr and L m ) and the mechanical parameters (H and F ) of a similar IM with similar R s and R r found in the literature have been used. The parameters used during the model and the sensitivity analyses are the following: R s = 0.196 Ohm R r = 0.0191 Ohm L ls = 0.0397 H L lr = 0.0397 H L m = 1.354 H H = 0.095 F = 0.0548 p = 1 (30) 112 Figure 5: The model state variables for a ±50% change of L m It can easily be seen that the IM model is nonlinear since there are products of two state variables in the following equations: • Eq.(27): �0 dr is multiplied by ! r • Eq.(28): �0 qr is multiplied by ! r • Eq.(29): �0 dr is multiplied by � qs • Eq.(29): �0 qr is multiplied by �0 ds For the local stability analysis we have to calculate the eigenvalues of the system. The examined equilibrium point is � qs = �0.744 Vs � ds = 1.0287 Vs �0 qr = �0.174 Vs (31) �0 dr = �0.087 Vs ! r = 48.477 s�1 The numerical value of the Jacobian of the nonlin- ear model (i.e. the state matrix of the locally linearized model) is as follows: 2 66664 �2.504 �50 2.4328 0 0 50 �2.504 0 2.4328 0 0.2370 0 �0.244 �1.522 0 0 0.2370 1.5222 �0.244 0 �8.617 17.077 0 0 �0.288 3 77775 . (32) The eigenvalues of the state matrix of the linearized sys- tems are: � 1,2 = �2.504 ± j49.98, � 3,4 = �0.243 ± j1.534, and � 5 = �0.288 (33) Figure 6: The model state variables for a ±50% change in R s It is apparent that the real parts of the eigenvalues are negative with small magnitudes. Parameter Sensitivity Analysis the sensitivity of the nonlinear model to the mutual induc- tance has been investigated. The steady state value of the system variables does not change (as is apparent in Fig.5) even for a considerably large change in L m . Sensitivity analysis of the inductances L lr , L ls of the stator and ro- tor, resistances of the stator R s and the damping constant F has also been investigated. In this investigation it can be seen that the values of the state variables were changed a bit, as shown in Fig.6. The analysis of the resistances of the rotor R0 r and the inertia H of the rotor showed that every value of the state variables changed significantly, as shown in Fig.7. As a final result of the sensitivity analysis, we can de- fine the following groups of parameters: • Not sensitive: Mutual inductance L m . Since the state space model of interest is insensitive in this respect, the values of this parameter cannot be reliably deter- mined from measurement data using any parameter estimation method. • Sensitive: These sensitive parameters are candidates for parameter estimation. – Less: inductances L lr and L ls , resistance of the stator R s and the damping constant F . – More: resistances of the rotor R0 r and the iner- tia H of the rotor. • Critically sensitive: none. 113 Figure 7: The model state variables for a ±50% change in H Conclusion The simple nonlinear dynamical model of an IM has been investigated in this paper. It has been shown that the model is locally asymptotically stable with regards to a physically meaningful equilibrium state. The effect of the controlled generator has been analyzed by simulation us- ing a traditional PI controller. It has been found that the controlled system is stable and can follow the set-point changes. Seven parameters of the model were selected for sensitivity analysis, and the sensitivity of the state vari- ables has been investigated. As a result, the parameters were partitioned into three groups. Based on the results presented here, a future aim is to estimate the parameters of the model for a real system from measurements. The sensitivity analysis enables us to select the candidates for estimation that are inductances L lr and L ls , resistances R s and R r , the damping constant F and the inertia H . An additional future aim is to develop a model for con- trol purposes and investigate different controllers before applying them on real systems. ACKNOWLEDGEMENT We acknowledge the financial support of this work by the Hungarian State and the European Union under the TAMOP-4.2.2.A-11/1/KONV-2012-0072 project. We ac- knowledge also the investigated IM of Grundfos Hungary Producer Ltd. SYMBOLS C 1 , C 2 , C 3 constant F damping constant � qs and � ds q and d components of the stator flux �0 qr and �0 dr q and d components of the reduced rotor flux H inertia constant i qs and i ds q and d components of the stator current i0 qr and i0 dr q and d components of the reduced rotor current L m mutual inductance p number of pole pairs R r and L r rotor resistance and inductance R0 r and L0 r reduced value of rotor resistance and inductance R s and L s stator resistance and inductance v d and v q q and d components of the stator voltage ! angular velocity of the magnetic field ! r angular velocity of the rotor T mech mechanical torque REFERENCES [1] VAS P.: Artifical-intelligence-Based Electrical Ma- chines and Drives, Oxford University Press, 1999 [2] VAS P.: Sensorless Vector and Direct Torque Con- trol. Oxford University Press, 1998. [3] ZHENG L., FLETCHER J.E., WILLIAMS B.W. HE X.: Dual-Plane Vector Control of a Five-Phase In- duction Machine for an Improved Flux Pattern, IEEE Transac. Ind. Elec., 2008, 55(5), 1996–2005 [4] LEVI E.: Impact of Iron Loss on Behavior of Vector Controlled Induction Machines, IEEE Transac. Ind. Appl., 1995, 31(6), 1287–1296 [5] HASEGAWA M., MATSUI K.: Robust Adaptive Full-Order Observer Design with Novel Adaptive Scheme for Speed Sensorless Vector Controlled In- duction Motors, Proc. 28th Annual Conf. Industrial Electronics Society, IEEE 2002, 2002, 1, 83–88 [6] GEYER T., PAPAFOTIOU G., MORARI M.: Model Predictive Direct Torque Control-Part I: Concept, Algorithm, and Analysis, IEEE Transac. Ind. Elec., 2009, 56(6), 1894–1905 [7] KRAUSE P.C., WASYNCZUK O., SUDHOFF S.D.: Analysis of Electric Machinery, IEEE Press, New York, 1995 [8] MOHAN N., UNDELAND T.M., ROBBINS W.P.: Power Electronics: Converters, Applications, and Design, John Wiley & Sons Inc., New York, 1995 [9] HATOS P., FODOR, A., MAGYAR A.: Parameter sensitivity analysis of an induction motor, Hung. J. Ind. Chem., 2011, 39(1), 157–161 Advertise upcoming meetings, conferences, workshops, new products or services; make public announcements, introduce your research laboratory on the pages of the Hungarian Journal of Industry and Chemistry Please contact us if interested! EDITORIAL OFFICE: UNIVERSITY OF PANNONIA P.O. BOX 158, VESZPRÉM H-8201 (HUNGARY) Tel.: +36 (88) 624-746, E-mail: hjic@almos.uni-pannon.hu; web: hjic.mk.uni-pannon.hu Felelős szerkesztő: Szilágyi Róbert Károly, PhD Kiadja: Pannon Egyetem, 8200 Veszprém, Egyetem u. 10. Levélcím: H-8201 Veszprém, Postafiók 158, Tel.: (88) 624-000 Felelős kiadó: a Pannon Egyetem, Mérnöki Kar dékánja