Parallel mathematical models of dynamic objects

The paper deals with the developing of the methodological backgrounds for the modeling and simulation of complex dynamical objects. Such backgrounds allow us to perform coordinate transformation and formulate the algorithm of its usage for transforming the serial mathematical model into parallel ones. This algorithm is based on partial fraction decomposition of the transfer function of a dynamic object. Usage of proposed algorithms is one of the ways to decrease calculation time and improve PC usage while a simulation is being performed. We prove our approach by considering the example of modeling and simulating of fourth order dynamical object with various eigenvalues. This example shows that developed parallel model is stable, well-convergent, and high-accuracy model. There is no defined any calculation errors between well-known serial model and proposed parallel one. Nevertheless, the proposed approach's usage allows us to reduce calculation time by more than 20% by using several CPU's cores while calculations are being performed.


Introduction
Nowadays different mathematical models are widely used in various scientific, economic, social, industrial and other applications [1]- [3].Various kinds of mathematical models make it possible such as to predict some events and actions [4] or make conclusions [5], and process some data and find interrelations between its various pieces [6]- [8].
One of the most important spheres of mathematical modeling and simulation is its scientific application [9]- [13].One can find that the benefits of mathematical models are very useful while novel objects, properties, and characteristics are being created [14]- [16] and/or it is impossible to perform some measurements.Moreover, tons of mathematical models allow studying objects which have not been created yet or other ways of their study are very dangerous and can cause accidents [17] immediately or after a while.
Class of dynamic models is a very wide class of mathematical models which allow us to study object's actions in process of time.Usage of such models makes it possible to study various dynamic and static characteristics and predict a behavior of the studied object.If one performs analysis of known dynamical models, he can find that these models are parallel if a considered object has parallel channels for transmitting energy from inputs to outputs, and these models are series otherwise.Many models of robotic devices [18] and some complex industrial processes are parallel models due to using many inputs and outputs but, nevertheless, most of the dynamic models are series ones.One of the most powerful tools for studying of object's dynamic is Simulink from MatLab software [19].This software is intended to solve various linear and nonlinear algebraic and differential equations which are given by state space

A R T I C L E I N F O A B S T R A C T
The paper deals with the developing of the methodological backgrounds for the modeling and simulation of complex dynamical objects.Such backgrounds allow us to perform coordinate transformation and formulate the algorithm of its usage for transforming the serial mathematical model into parallel ones.This algorithm is based on partial fraction decomposition of the transfer function of a dynamic object.Usage of proposed algorithms is one of the ways to decrease calculation time and improve PC usage while a simulation is being performed.We prove our approach by considering the example of modeling and simulating of fourth order dynamical object with various eigenvalues.This example shows that developed parallel model is stable, well-convergent, and high-accuracy model.There is no defined any calculation errors between well-known serial model and proposed parallel one.Nevertheless, the proposed approach's usage allows us to reduce calculation time by more than 20% by using several CPU's cores while calculations are being performed.equations as well as transfer functions.One can conclude that numerical solutions of differential equations are performed in a serial way.This fact determines the main drawback of differential equations' solvers which are integrated into Simulink.These solvers define solution of differential equations one after the other and use for this action only one core from multi-core CPU.This drawback increases calculation time dramatically and causes under-utilization of modern PC while a dynamic of the object is being studied.One can find similar problems with other mathematical software as well.
We offer to decrease calculation time and improve CPU usage by developing a scientific background for developing the parallel model that can be simulated by using different cores of CPU, CPUs, and PCs.Parallel computations are widely used while scientific research is being carried out [20], [21] but no parallel dynamic models were found while scientific publications were analyzed.Well-know Parallel Computational Toolbox which is part of Matlab does not allow to decompose model into small parallel pieces and calculate them in a parallel way.

Transformation in Continuous Time Domain
We consider dynamic of a generalized dynamical object given by the following linear state space equations in matrix form [22].

  
where is the -sized state space vector, is the -sized output vector, is the -sized vector of input efforts, A is the -sized state matrix, B is the -sized input matrix, C is the -sized output matrix and D is the -sized feed through matrix.
It is easy to define matrix transfer function by using (1) Here, E is an identity matrix.
If one takes into account the matrix inversion rule, (2) can be rewritten as follows where is the characteristic polynomial of matrix A, is the adjugate matrix to matrix    here is the (i,j) minor of matrix.
Finally, we reduce (3) to a common denominator and write down following matrix transfer function Analysis of (5) makes it possible to formulate the following statements: Statement 1. Generalized matrix transfer function of a considered linear object consists of two parts.The first one defines the controlled dynamic of the considered object and the second one defines the influence of control signal on observed state variables.

Statement 2.
If state space equations have D matrix with non-zero elements, the order of numerator's and denominator's polynomials equal to object order.We assume that A, B, C, D matrices are well-defined and a result of numerator's (5) calculation can be rewritten by using matrix polynomial where is the some -sized matrix.
The denominator of ( 5) can be given in a similar way Expressions ( 6) and (7) give us possibility to rewrite (5)

  
Transfer function (8) makes it possible to write down following differential equation which we call as matrix differential equation in canonical form.

  
It is quite clear that equation (10) can have zero, real and complex eigenvalues.Moreover, these eigenvalues can be k-th multiple eigenvalues.Now we assume that (10) has k0 multiple zero eigenvalues, k11 real eigenvalues of order k12 multiplicity and k21 complex conjugate eigenvalues of order k12 multiplicity.This assumption allows us to rewrite denominator (8) as follows where is the i-th real eigenvalue, and are the real and imaginary parts of j-th complex conjugate eigenvalue, I is the imaginary unit.
Expression (11) allows us to rewrite (8) as follows   and decompose second multiplier in following way [23]   where are some coefficients.
One can find these coefficients by solving the following equation We offer to use (13) for rewriting (12)   We call (15) as parallel matrix transfer function.
Now we assume the existence of some virtual state variable and define state vector in the following way

  
This assumption makes it possible to define (15) as the sum of subsystem transfer functions It is clearly understood that the order of each summand (17) defines its order Expressions ( 15) -( 19) allows us to make the following statements: Statement 3. Usage of partial fraction decomposition allows us to rewrite object's transfer function into parallel way by defining virtual parallel channels.Order each of them less than the order of whole transfer function.
Statement 4. One can study the dynamic of the generalized linear object in an easy way by studying the dynamic of all parallel channels.
We suggest determination of state space equations for a virtual parallelized object by using ( 16)-( 18) as follows

 
We call (30) as parallel state space equations.
It is clearly understood that (30) contains n independent differential equations which can be solved separately analytically or numerically.If one solves them in a numerical way he can use different cores and CPUs or even different PCs to solve (30).When values of virtual state space vectors are obtained, one can store and sum them for getting values of vector.This is the main benefit of the proposed approach which can be implemented by means of computer technology in an easy way.

State space equations in the discrete time domain
There are tons of various numerical methods for solving ordinary differential equations which differs with various approximations of the derivative operator.The common feature of these methods is a usage of state variables in previous, current, and future time.
Thus, we consider the generalized approximation of derivative operator where is a shift operator, is a number of previous time's samples, is a number of future time's samples.

Usage of (21) allows us to rewrite (20) in discrete time domain
Equations ( 22) allow us to formulate the following statement.
Statement 5. Solver, which calculates dynamic for every parallel channel, operates with pre and post information about virtual state variables' values and control efforts in the general case.One can use (22) for defining virtual and real state variables for various sample time T.

Parallel model of DC electric drive in continuous time
It is well-known fact that one can describe dynamic of DC electric drive by the following differential equations which are given by using the relative units

  
where coefficients are defined as follows here J is a rotor inertia, R is an armature resistance, c is a back-emf constant, L is an armature inductance, is power converter coefficient.
We rewrite (23) into matrix form One can find that this condition is true for the most industrial application, so we assume that (28) has two complex conjugate eigenvalues

  
This fact makes it possible to perform for (28) partial fraction decomposition We offer to perform for (33) some algebraic transformations and simplifications which allows us to write down equations for defining unknown coefficients and define them in the following way   We use (33) to write down parallel state space equations as follows It is clearly understood that the transformed mathematical model has three parallel channels.Two of them have first order and the third one is a dynamical object of second order.In such a way we replace one complex fourth order dynamical object with three simple dynamical objects.

Parallel model DC electric drive in the discrete time domain
Now we transform (37) into the discrete time domain by applying bilinear transformation which allows us to approximate derivative operator in the following way and rewrite (37) One can use (27) or similar approximations of the derivative operator for obtaining the mathematical model of DC electric drive in the discrete time domain.Right-hand expressions of (40) always depend on previous values of virtual state variables, but this dependence is different for various approximations.
It is well-known fact that these approximations define stability and accuracy of numerical solution (37).From a programming viewpoint, it is easy to solve so-called mesh equations.

Parallel mesh model of DC electric drive
Now we rewrite (40) as mesh equations One can use (43) and (45) for defining output variable of considered DC electric drive.

Simulation's Results
Results of numerical solutions ( 23) and (37) are shown in the Fig. 1.
We obtain these results in the discrete time domain with sample time and parameters .Simulation time for the serial model is 0.4299s while one CPU core is being used.First equation of the parallel model is solved in 0.0878s, the second equation is solved in 0.1347s, and the last two equations are solved in 0.2615s.It needs 0.0651s for summing the results of each parallel channel.So, simulation time for the parallel model is 0.3266s by using 3 CPU cores.Thus, we decrease simulation time by 24% while we use parallel model and its multi-core program implementation.Analysis of simulation's results makes it possible to conclude that serial and parallel models allow us to get the same results.Moreover, numerical solutions for both models coincide with great accuracy which is less than sample time.

A Generalized Algorithm for Transformation of Serial Models into Parallel Ones
Above-given theoretical and numerical results allow us to formulate an algorithm for transformation serial mathematical continuous time model into parallel discrete time one.There are two operations which can be concluded:  State space equations for new virtual state variables are written down.
2) Operations which are performed in the discrete time domain, which is containing:  Approximation formula for the derivative operator is chosen.
 State space equations are rewritten into the discrete time domain.
 Mesh equations are written down.

Conclusion
The mathematical model of the generalized linear dynamical system can be transformed into parallel form by performing partial fraction decomposition for its transfer function.Usage of this transformation allows us to simplify the order of considered dynamical system and study its dynamical and static ones if the following condition is true  

Fig. 1 .
Fig. 1. Results of simulations serial and parallel models

1 )
Operations which are performed in the continuous time domain, which has produced:  State space equations for the considered object are written down. The transfer function is defined. Characteristic polynomial is analyzed and eigenvalues are defined. Partial fraction decomposition for transfer function is performed.