Nonlinear Continuous System Identification by Means of Multiple Integration II

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.


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) and an output equation (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 followed by an algebraic block ( ) u, .
2. There exists an inverse function ( ) g y y u, of the output function ( ) u, 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 ( ) u, 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 t where The functions f i are supposed to be linear in parameters a i, j , i.e., to have the form where g i, j are known (generally nonlinear) functions of the measured data and a i, j (generally unknown) constants. Example:

Identification algorithm
For the sake of simplicity, let us denote a multiple integral of a time-dependent variable x(t) as ( ) , , 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 x 1 (t) in Fig. 1 can be expressed as From (3), ( 4), (7) follows Let us denote the weighted sum of multiple determined integrals (6 ) (after enumeration of this expression we obtain a scalar constant) as From ( 5), ( 11), ( 12) follows the linear algebraic equation for the unknown coefficients a i, j The weighted sums C of determined integrals (12) of the (generally nonlinear) functions ( ) , u represent in (13) the known linear coefficients of the algebraic linear equation for the unknown constants a i,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 a i,j , for example by the least squares method.To obtain a unique solution, one or more constants a i,j must be known.These constants are then put on the right side of system (13).The MATLAB program MI, described in part 4, was designed to enable such computation.

Examples
A. A linear system, described by a transfer function 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 & ; , , , & and the output equation is ( ) Linear algebraic equations for calculating the unknown system parameters are An example of an identification algorithm for a linear system 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 w.The output is dynamo voltage e. (The dynamo is supposed without any external load).The differential equation for the system dynamics is ( ) where e = w F is the induced voltage, F is the machine flux, The inverse magnetization curve i e (F) is supposed in the form Let us define the variables and constants to correspond to the symbols introduced in part 1.The negative inverse output function is then ( ) The algebraic equation (13) for computing the unknown parameters is ( ) In original notation:    (24)

Identification program MI and its use
The identification program (available in [4]) is written for one-shot identification (identification off-line).The basic condition 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 necessity 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 future model.(These frequencies often lie around the critical frequency w 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]) T min £ p/w 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 T min .To obtain one algebraic equation ( 13), a time sequence (n+2)T of data is necessary.In the program, the time origin of the data for a new equation ( 13) is shifted by T/2.This implies that from the time period t max , approximately n e (T) = 2 ( t max -(n+2) T) /T+1= 2t max / 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 n e (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 further details).It is recommended to choose the integration periods in a geometric series beginning by T min = p/w 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, consequently, the error of the calculated parameters.Let us suppose in the following that the model contains m independent nonlinear blocks g i,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 integrated 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 beginning of the sampling period; data is a matrix consisting of m column vectors of (generally) nonlinear functions of measured data, corresponding 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 members in the columns correspond to the index (i + 1), i.e., to the number of integrators between the block g i,j and the final state equation x 1 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, specifically according to the order of its nonzero elements.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 + + + 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 identification, column data vectors of measured data u(t) and y(t) will be necessary.Parameter a 0 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 x 1 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

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].

)Fig. 1 :
Fig. 1: State space diagram of the identified model

F
is the voltage induced by the change of the machine flux in the excitation winding, and Ri e is the voltage drop on the resistance in the excitation circuit.

Fig. 2 :
Fig. 2: State-space diagram of the linear system

Fig. 4 .
Fig. 4. State-space diagram of the parallel dynamo , 1) = -1.If the structure matrix has the form c = the result (output vector) parameters = [a 1 , a 2 , b 1 , b 2 ].Parameter a 1 corresponds to the member -2 in c.Consequently, a 2 corresponds to -3, b 1 to 2 and b 2 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 will be parameters = [b 1 , a 0 , a 1 , a 2 ].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 u 1 , u 2 and y (resistance in the excitation circuit R, angular velocity w, output voltage of the dynamo e), and if we want to obtain the program output vector parameters = [a 1, 1 , a 1, 2 , a 1, 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./w,e, R.*e./w, R.*((e./w).Ù 3)] and the output vector of results corresponds to the expression parameters = [1/N, a/N, b/N], see (20).