AP1_01.vp 64 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 41 No.1/2001 1 Basic definitions Let us take a nonlinear continuous time-invariant SISO (Single Input – Single Output) system, described by the state equations (see also Fig.1) � � � � � � � �� � � � � � � �� � � , ; , , , � , ; x x x i i+ n t t f t y t i n t f t y t i n i n � � � � � � 1 1 2 1u u � (1) and an output equation � � � �� � � � � �� �y t g g t x t ty u� ��1 1u u, . (2) (The state-space representation is substantially changed in comparison with [1].) Let us suppose that u and y are the only measurable quantities in the system. The state-space diagram must fulfill the following rules: 1. The basic string of the scheme consists of n integrators fol- lowed by an algebraic block � �g xy �1 1u, . 2. There exists an inverse function � �g yy u, of the output function � �g xy �1 1u, from Eq. ( 2 ). 3. All the direct paths of the scheme lead from the output source u on the basic string. 4. All the feedback loops of the scheme lead via the output algebraic block � �g xy �1 1u, and none of them is algebraic, i.e., all of them contain at least one integrator. Then we can defined � � � � � �� �x t f t y t1 0� � u , (3) where � � � �� � � �� � � � � �� �f t y t g t g t y tu y0 u u u, ,� � . (4) The functions fi are supposed to be linear in parameters ai, j , i.e., to have the form � � � �� � � � � �� �f t y t a g t y t i ni i j i j j mi u u, , ; , , ,, ,� � � 1 0 1 � , (5) where gi, j are known (generally nonlinear) functions of the measured data and ai, j (generally unknown) constants. Example: � � � �f y a u y a y a u u y3 3 1 1 3 2 3 3 3 1 2u, cos, , ,� � � � � . 2 Identification algorithm For the sake of simplicity, let us denote a multiple integral of a time-dependent variable x(t) as � � � �� � x i t x t t t t t i x tt tt i i i , , ; , , � � � �00 10 1 20 1 13 0 0 d d d d� � � � �t x t i� � 0. (6) Note that (6) corresponds to a time response of a chain of i integrators in series, with zero initial conditions, excited by the input function x(t). The first determined integral from 0 to t of the last state-space variable x1(t) in Fig.1 can be expressed as � � � � � � � �� �� �� �x t x t i f t y t i ti i i i n , , ! , , ,1 0 1 1 � � � � � � � � � � � � u . (7) From (3), (4), (7) follows � � � � � �� �� �x t i f t y t i ti i i n i i n 0 1 0 1 0 ! , , , � � � � �� u . (8) Now let us integrate (8) for kT, k = 1, 2, …, n + 2 and mul- tiply both sides of the k-th equation by a binomial coefficient � �� � � � �� � � ��1 2 k k n . Let us make a sum of the resulting (n + 2) equations to obtain Nonlinear Continuous System Identification by Means of Multiple Integration II J. John This paper presents a new modification of the multiple integration method [1, 2, 3] for continuous nonlinear SISO system identification from measured input – output data. The model structure is changed compared with [1]. This change enables more sophisticated systems to be identified. The resulting MATLAB program is available in [4]. As was stated in [1], there is no need to reach a steady state of the identified system. The algorithm also automatically filters the measured data with respect to low frequency drifts and offsets, and offers the user a potent tool for selecting the frequency range of validity of the obtained model. Keywords: continuous system identification, multiple integration. u y Fig. 1: State space diagram of the identified model � � � � � � � � � �� � � � � � �� �� � � � � � � � 1 2 0 1 1 2 1 k k n i i i n i k n x kT i f t y t i ! , ,� u� �, .t i n � � � � � � � � � 0 0 (9) Because � �� � � �� �� � � 1 0 1 2 1 1 k k n ik m k i m; , , ,� , (10) the generally non-measurable initial conditions xi(0) disap- pear from (9) and we obtain � � � � � �� �� �� � � � �� �� � � � � � � � � � � �� 1 2 1 01 2 k i i n k n k n f t y t i t� u , , , 0 . (11) Let us denote the weighted sum of multiple determined integrals (6 ) (after enumeration of this expression we obtain a scalar constant) as � � � � � �� �x i n T k n x i kTk k n , , , , , � � � � �� �� � � � 1 2 1 1 2 . (12) From (5), (11), (12) follows the linear algebraic equation for the unknown coefficients ai, j � � � �� �� �a g t y t i n Ti j i j j m i= n i , , , , , ,� u �� 10 0 (13) The weighted sums C of determined integrals (12) of the (generally nonlinear) functions � �g yi j, ,u represent in (13) the known linear coefficients of the algebraic linear equation for the unknown constants ai,j. By shifting the time origin of the input – output data and repeating algorithm (13) we obtain another algebraic equation. After repeating this process enough times we obtain a system of algebraic equations that can be solved for the unknown constants ai,j, for example by the least squares method. To obtain a unique solution, one or more constants ai, j must be known. These constants are then put on the right side of system (13). The MATLAB pro- gram MI, described in part 4, was designed to enable such computation. 3 Examples A. A linear system, described by a transfer function � � � � Y s U s b s b s b s b s b a s a s n n i n i n n n n � � � � � � � � � � � � � 0 1 1 1 0 1 1 � � � a s a s ai n i n n � �� � �� 1 (14) can be described by the state-space diagram in Fig.2. From the canonical forms this is the unique scheme corresponding to the structure shown in Fig. 1. State-space equations of the system are � ; , , , � x x b u a y i n x b u a y i i i i n n n � � � � �1 1 2 1� (15) and the output equation is � �y x b u a � � �1 0 0 1 . (16) Linear algebraic equations for calculating the unknown system parameters are � � � �� �b u i n T a y i n Ti i i n � �, , , , , ,� � 0 0 . (17) An example of an identification algorithm for a linear sys- tem of the fourth order can be found in [1]. B. A parallel dynamo (Fig. 3) is a nonlinear system with two input variables: resistance of the excitation winding R and angular velocity of rotation �. The output is dynamo volt- age e. (The dynamo is supposed without any external load). The differential equation for the system dynamics is � �e =Ri N t e �� � � �d d , (18) where e � � � is the induced voltage, � is the machine flux, N t d d � is the voltage induced by the change of the machine flux in the excitation winding, and Rie is the voltage drop on the resistance in the excitation circuit. The inverse magnetization curve ie(�) is supposed in the form � �i a be � � � � � . (19) Let us define the variables and constants to correspond to the symbols introduced in part 1. u R u x y e u x a N a a N a b N 1 2 1 2 1 1 1 1 2 1 3 1 � ; ; ; ; ; ; ., , , � � (20) The state-space diagram in Fig. 4 corresponds to (18). It does not fulfill the condition of closing the feedback loops via © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 65 Acta Polytechnica Vol. 41 No.1/2001 u y Fig. 2: State-space diagram of the linear system Fig. 3: Parallel dynamo the measurable output variable y (item 4, page 64). The dia- gram must therefore be restructured, as in the Fig. 5. Then it corresponds to a state-space equation � , , ,x a y a u y u a u y u 1 1 1 1 2 1 2 1 3 1 2 3 � � � � � �� � � (21) The negative inverse output function is then �f u y y u 0 2 , � � . (22) The algebraic equation (13) for computing the unknown parameters is �� � �� � � �� � � � � � y u T a y T a u y u T 2 1 1 1 2 1 2 0 1 1 1 1 1, , , , , , , , ,, , � � �� � � � � � � � �� � � � � � � � � �a u y u T1 3 1 2 3 1 1 0, , , ,� (23) In original notation: � � � � � � � � � � � � � � � � � � � � � � � 2 1 2 0 0 2 1 1 0 2 e t t t e t t t N e t t T T t � � d d d d � � � � t e t t t a N R t e t t t tTT 2 1 1 0 2 0 2 0 1 1 1 2 2 � � � � � � � � � � � � � ��� d d d � � � � 1 0 2 1 1 1 1 0 2 0 2 0 2 2t tTT t R t e t t t t b N � ��� � � � � � � � � � � � � � d d d � � � � � � � 2 1 1 1 1 0 2 1 1 1 2 R t e t t t t R t e t t t t � � � � � � � � � � � � � � � � 3 3 d d d 1 0 2 0 2 0 2 0 tTT t��� � � � � � � � � � �d . (24) 4 Identification program MI and its use The identification program (available in [4]) is written for one-shot identification (identification off-line). The basic con- dition for using the program is to find a proper structure of the model, corresponding to the state-space diagram in Fig. 1. This is not always a simple task and sometimes we need to transform the physical reality, as we did it in the case of the parallel dynamo in item B of the previous part. The other ne- cessity is the existence of relevant data u, y in the form of time vectors. The data must contain a considerable content of frequencies interesting from the point of view of the fu- ture model. (These frequencies often lie around the critical frequency �k, i.e., during the measurement, at least for a short time, the identified system should oscillate in phase opposition.) The shortest integration period T, which is one of the items of the program input data, has to be (see [1]) Tmin � �/�k. By Shannon’s theorem, the data ought to have at least ten samples in the critical period. This implies for the sample period h � 0.2 Tmin. To obtain one algebraic equation (13), a time sequence (n +2)T of data is necessary. In the pro- gram, the time origin of the data for a new equation (13) is shifted by T/2. This implies that from the time period tmax, ap- proximately ne(T) = 2 ( tmax – (n +2 ) T) /T+1= 2tmax / T – 2 n –3 equations will be obtained (an approximation is given by rounding). More than one integration period can be put for one program call, and for each of them its own system of ne(T) algebraic equations (13) will be obtained. Before the final computing of the parameters, all these systems are joined together and the parameters are calculated from this joint system by the method of least squares. By a proper choice of the integration period vector T we can freely select the filtration properties of the kernels (see [1] or [4] for fur- ther details). It is recommended to choose the integration periods in a geometric series beginning by Tmin = �/�k and with quotient q � 1.5�4. The program displays on the screen the numbers of equations used for the given integration period T, the sum of these numbers, the root mean square error of the complete equation system, and the singular numbers of its matrix. The singular numbers can serve as a measure of the system conditionality number and, con- sequently, the error of the calculated parameters. Let us suppose in the following that the model contains m independ- ent nonlinear blocks gi,j from the equation (5). These blocks must be programmed as vector functions of the measured data vectors (by the period convention of MATLAB). The program is called by a statement parameters = mi(h,T,staircase,data,c) where h is the sampling periods of the measured data T is the row vector of the sampling periods staircase is a row m-vector of data indicators; for starcase(i) = 0 the corresponding data will be in- tegrated as continuous, by Euler’s method, for staircase(i) = 1 the corresponding data will be integrated as stepwise with a valid value at the end of the sampling period, for staircase(i) = �1 the corresponding data will be integrated as stepwise with a valid value at the be- ginning of the sampling period; 66 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 41 No.1/2001 Fig. 4. State-space diagram of the parallel dynamo Fig. 5: Restructured state-space diagram of the parallel dynamo data is a matrix consisting of m column vectors of (gen- erally) nonlinear functions of measured data, cor- responding to the functions from equation (5); c is a matrix describing the model structure. It consists of m columns, corresponding to the columns of the matrix data. Absolute values of the nonzero mem- bers in the columns correspond to the index (i + 1), i.e., to the number of integrators between the block gi,j and the final state equation x1 at the end of the integrator chain in Fig. 1, augmented by one. The sign of the member corresponds to the sign of the connection. The coefficient c (1, 1) corresponds obligatorily to the known parameter (here one), and the corresponding sum of weighted integrals will be put with inverted sign on the right side of equation (13). We must complete matrix c by zero members to a rectangular shape. These members are of no significance for the computation; parameters is the output row vector. It is ordered in accordance with the columns of matrix c, specifi- cally according to the order of its nonzero ele- ments. Unitary coefficient corresponding to c(1, 1) does not appear in the output vector. Examples A. The linear system with a transfer function b s b s a s a 1 2 2 1 2 � � � corresponds to the state-space diagram in Fig. 2 for n = 2. The corresponding algebraic equations (13), built from the weighted integrals, will have the form of (17). For identifica- tion, column data vectors of measured data u(t) and y(t) will be necessary. Parameter a0 is selected as unitary in our transfer function. In the state-space diagram, Fig. 2, this parameter is connected with the last state-space variable x1 only by an algebraic connection (without integrators), and in the equations it appears with negative sign. Due to this reasons the matrix data will be ordered as [y, u], and in the structure matrix c the first member will be c (1, 1) = �1. If the structure matrix has the form c � � � � � � � � � � � � � � 1 2 2 3 3 0 , we obtain the result (output vector) parameters = [a1, a2, b1, b2]. Parameter a1 cor- responds to the member –2 in c. Consequently, a2 cor- responds to –3, b1 to 2 and b2 to 3. Parameter 1 by s 2 in the transfer function denominator corresponds to the member –1 in c. For the same measured data we can obtain the resulting transfer function in the form b s a s a s a 1 0 2 1 2 1� � � with the inputs data = [u, y], c � � � � � � � � � � � � � � 3 1 2 2 0 3 and the result will be parameters = [b1, a0, a1, a2]. B. Parallel dynamo – see Fig. 3, 4, 5 and Eqs. (19–25). If we have at our disposal the time vectors of the measured data u1, u2 and y (resistance in the excitation circuit R, angular velocity �, output voltage of the dynamo e), and if we want to obtain the program output vector parameters = [a1, 1, a1, 2, a1, 3], the input data matrix must be (description in MATLAB) data=[� y./u2, y, u1.*y./u2, u1.*((y./u2).�3)] and the matrix of description of the model structure must be c=[�1, 2, �2, �2]. In the original notation the input data matrix is data=[e./�, e, R.*e./�, R.*((e./�).�3)] and the output vector of results corresponds to the expression parameters = [1/N, a/N, b/N], see (20). 5 Conclusions A new modification of the methods presented in [1, 2, 3] has been presented, including the corresponding program. In comparison to [1, 2, 3], more general nonlinearities are permitted. A complete description of the method, more examples and corresponding programs can be found in [4]. References [1] John, J., Rauch, V.: Nonlinear continuous system identifica- tion by means of multiple integration. Acta Polytechnica, Vol. 39 (1999), No. 1, pp. 55–62 [2] Maletínský, V.: I-i-p Identifikation kontinuierlicher dyna- mischer Prozesse. Abhandlung zur Erlangung des Titels eines Doktors der Technischen Wissenschaften der ETH Zürich, (In German: I-i-p Identification of Continuous Dy- namic Systems. PhD thesis, ETH Zürich.), SSS Zürich 1978 [3] Skočdopole, J.: Nelineární vícenásobná integrace. Kandi- dátská disertační práce FEL ČVUT 1981, (In Czech: Nonlinear multiple integration.) PhD thesis, FEE CTU Prague 1981 [4] See: http://dce.felk.cvut.cz/sri2/index.htm Doc. Ing. Jan John, CSc. Department of Control Engineering phone: +420 2 2435 7597 e-mail: john@control.felk.cvut.cz Czech Technical University in Prague Faculty of Electrical Engineering Technická 2, 166 27 Praha 6, Czech Republic © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 67 Acta Polytechnica Vol. 41 No.1/2001 Table of Contents Thermal and Hygric Expansion of High Performance Concrete 2 J. Toman, R. Èerný Temperature and Moisture Dependence of the Specific Heat of High Performance Concrete 5 J. Toman, R. Èerný Thermal Conductivity of High Performance Concrete in Wide Temperature and Moisture Ranges 8 J. Toman, R. Èerný Faraday Cup for Electron Flux Measurements on the Microtron MT 25 11 M. Vognar, È. Šimánì, D. Chvátil Some Aspects of Profiling Electron Fields for Irradiation by Scattering on Foils 14 È. Šimánì, M. Vognar, D. Chvátil Critical Length of a Column in View of Bifurcation Theory 18 M. Kopáèková Erosion Wear of Axial Flow Impellers in a Solid-liquid Suspension 23 I. Foøt, J. Medek, F. Ambros A Small Transfer and Distribution System for Liquid Nitrogen 29 È. Šimánì, M. Vognar, J. Køíž, V. Nìmec Lewatit S100 in Drinking Water Treatment for Ammonia Removal 31 H. M. Abd El-Hady, A. Grünwald, K. Vlèková, J. Zeithammerová Evaluation of Water Resistance and Diffusion Properties of Paint Materials 34 J. Drchalová, J. Podìbradská, J. Madìra, R. Èerný Management of People by Managers with a Technical Background – Research Results 38 D. Dobrovská Clinoptilolite in Drinking Water Treatment for Ammonia Removal 41 H. M. Abd El-Hady, A. Grünwald, K. Vlèková, J. Zeithammerová Application of Morphological Analysis Methodology in Architectural Design 46 A. Prokopska Experimental Investigations on Drying of Porous Media Using Infrared Radiation 55 A. K. Haghi Dynamic Effect of Discharge Flow of a Rushton Turbine Impeller on a Radial Baffle 58 J. Kratìna, I. Foøt, O. Brùha Nonlinear Continuous System Identification by Means of Multiple Integration II 64 J. John