Abstract-Currently, devices of synchronized vector measurements (Wide-Area Monitoring Systems-WAMS) are becoming more widespread in power systems. With a high accuracy, they measure complex current and voltage values and also frequency in places of their installation. Obtained data allow to create an objective picture of power system state. This paper shows the possibility of using WAMS measurements for parameters estimation of synchronous machines (SM) and their automatic excitation controllers (AEC). At present for a variety of tasks widely used simulation modeling of synchronous machines (SM) in software and hardware-software packages. In this case there are often difficulties in forming accurate models of real SM due to the lack of raw data. In view with the foregoing, the urgent task is to create a tool for clarifying and determining the SM model parameters. This tool can be formed based on various optimization algorithms (including genetic algorithm). For definition of SM parameters sharing WAMS devices and optimization algorithms is perspective and actual research direction. To solve this issue, can be used wellknown differential equations system of Park-Gorev. All phases currents and stator voltage, voltage and current of excitation, network frequency are easily detected through WAMS and are known values. The availability of data on SM electric values combined with using differential equations system of Park-Gorev allows to determine of SM parameters with prescribed periodicity. As a result of further development of the proposed technique of determining SM parameters it can be the basis for SM parameters monitoring system according to the data obtained from WAMS in real time. Another relevant application of the developed technique is the solution of the issue of power system equivalenting or its part according to the data obtained from WAMS. Working of synchronous machines is impossible without automatic excitation control. This paper shows the possibility of using WAMS measurements for estimate AEC parameters of synchronous generator for different AEC types from various manufacturers in ideal and real conditions. By using the simulation in the hardware and software system RTDS (Real-Time Digital Simulator), WAMS measurements have been obtained for different modes of AEC operation. The developed technique is based on the processing of obtained measurements by different optimization algorithm. As a result of the comparative analysis, the most suitable type of # I. INTRODUCTION Currently, phasor measurement units (PMU) are becoming more widespread in power systems. With a high accuracy, they measure complex current and voltage values in the installation sites of PMU. Obtained data allow creating an objective picture of power system state. t this time, in carrying out various researches, simulation of synchronous generators (SG) is widely used. The possibilities and, consequently, the list of parameters that need to be set in the SG block model, in a variety of Program Complexes (PCs) and Program Apparatus Complexes (PACs) are different. At that, a situation when there are no values of various parameters required for setting in the SG models often happens. This factor has a negative influence in forming the model and its further use. To determine the SG model parameters both classical methods described in [1] and different from them can be used. For example, different optimization algorithms can be successfully used for this, that will increase the accuracy of simulation SG. an object of simulation can be a real SG or SG model, made in another PC, PAC. Necessary data to determine the parameters of the real SG can be obtained by PMU usage. To determine the Xd parameter by the traditional method in accordance with [1,3] is required to construct open-circuit characteristic (OCC) and three-phase shortcircuit characteristic (SCC). These characteristics must be built on the same plane, and then for a certain value of the excitation current (is taken equal to 1,0) Xd parameter as the ratio is a voltage (for OCC) to a current (in SCC). The resulting parameters value will be saturated. You must perform to obtain an unsaturated Xd value the same steps using straightening unsaturated OCC. OCC and SCC, built for the SG model, constructed in PC Matlab, shown in Fig 1. Xd definition by using the optimization algorithm implies the comparison and the minimization of the difference between the OCC, belonging to the object of simulation, and OCC obtained in formed model, by automatically selecting the Xd parameter value. At that time, the OCC, belonging to the synchronous generator model, is benchmark characterstic, and OCC, belonging to the newly formed SG model, are automatically reset for each new set Xd value, as long as it does not match the benchmark OCC with a given accuracy. To determine the X'd, X''d, T'd, T''d parameters by the traditional method in accordance with [1] is required to hold the experience sudden three-phase short circuit (TSC). Fixed open circuit voltage immediately prior the short circuit (SC) and the phase currents of the stator SG. X'd is defined as the ratio of the open circuit voltage to the initial value of the periodic component of the short-circuit current net subtransient component. X''d is defined as the ratio of the open circuit voltage to the initial value of the periodic component of the short-circuit current. T'd is defined as the time during which the transient component of the stator current decreases to its initial value of 0.368. T''d is defined as the time during which the subtransient stator current decreases to its initial value of 0.368. To determine the X'd, X''d, T'd, T''d parameters by using optimization algorithms an optimization function, whose arguments are all required to parameters, is formed. On the SG model the experience sudden TSC is held and an oscillogram of phase stator current is fixed. From fixed oscillograms stands out periodic component of the current is defined as the halfalgebraic ordinates of the upper and lower envelopes amplitude (shown in Fig. 2). After that the difference between the resulting curve and standart is minimized. III. DETERMINATION OF RS, TA, X 2 The Rs and Ta parameters should be determined from the experience of a sudden three-phase short circuit on the attenuation of the periodic component of the current in the excitation circuit, by the traditional method in accordance with [1]. For this the oscillogram of the excitation current is fixed at the shorted SG phase stator winding, which works in an idling. From the fixed oscillogram, the periodic current component stands out and is plotted in semi-logarithmic coordinates. Value of the time constant of the aperiodic stator current component Ta is defined as the time during which a periodic component of the current drive circuit decreased to its initial value of 0.368. Rs = 0.8?X 2 / (? c ?T a ) (1) Where X 2 is an unsaturated parameter and can be determined by the equation [2]: X 2 = 1.1?X''d(2) The approach to the determination of the Rs and Ta parameters by using optimization algorithms is similar to the one, which was used to determine the X'd, T'd parameters. Optimization function, whose argument is the desired parameter Rs, is formed. On the SG model the experience sudden TSC is held. From a fixed excitation current oscillogram the periodic component stands out. After that the difference between the resulting curve and benchmark is minimized. The subsequent determination of the Ta time constant is produced by the equation (1). At the same time X 2 is defined by ( 2), wherein the subtransient inductive reactance X''d has already been defined by the optimization method. It should be noted, that all data about necessary electrical quantities can be obtained by PMU usage. # IV. DETERMINATION OF t'd0 For the determination of the transient time constant along a longitudinal d-axis at the open stator winding by the traditional method in accordance with [1], a voltage recovery method can be used. The difference between the steady-state value voltage and recovery voltage U ( -U, defined by their amplitudes, is plotted in semi-logarithmic coordinates. Determination of the T'd0 parameter by the voltage recovery method consists in determining time during which a transient voltage component Î?"U' decreases to its initial value of 0.368. The procedure for determining T'd0 parameter by the optimization algorithms similar to the procedure, which was performed in determining Xd parameter. At the same time is held the voltage recovery experience on the SG model and the stator voltage recovering oscillogram is also fixed, then from fixed oscillogram stands out the envelope and the difference between the benchmark curve and the obtained on the SG model is minimized. # Global # DETERMINATION OF SG STATOR WINDING LEAKAGE INDUCTIVE REACTANCE XL XL parameter according to the [1] can be determined graphically by OCC, SCC and the point of load characteristics, which corresponds to the nominal values of voltage and current in the stator overexcitation mode (shown in Fig. 3). # Figure 3: Graphical determination of the SG stator winding leakage inductive reactance Xl To determine the stator winding resistance, the equation can be used [2]: ?) To the left of the point A parallel to axis of abscissas AF segment is laid, equal to the value of the excitation current to be determined at the nominal stator current on the steady TSC characteristic. From point F the line parallel to the initial part of the OCC to intersect the latter at a point H is laid. Perpendicular from the point H on the line AF there is a voltage drop across the leakage inductive reactance Xl at the nominal stator current. Calculated inductive reactance Xl is determined by the equation: Xl = U HG (6) Plotted on the Fig. 3 characteristics can be obtained using PMU data (about changes of electrical quantities). Determinations of the Xl parameter by the optimization method is performed similarly to the synchronous inductive reactance along a longitudinal d-axis Xd. The only difference is that the argument of the optimization function is the stator winding leakage inductive reactance resistance Xl. # VI. DETERMINATION OF tj, h The nominal time of the machine acceleration Tj and the storage-energy constant H can be determined by the run-out idling method, by the traditional method according to [1]. Rate speed test machine is set above the nominal, then the power source is disconnected. At the same time the SG rate speed change is fixed, then the time ?t, during which the machine changes rate speed in the range of Î?"n, is determined. To determine the H, Tj parameters by the optimization method the SG rate speed changes oscillogram, which is fixed during the run-out idling experiment, should be used. The optimization function, whose argument is a storage-energy constant H, is formed. On the model of SG the run-out idling experiment is held, and then the difference between the resulting curve and standard is minimized. Tj parameters can be determined by the equation, wherein the storage-energy constant H has already been defined by the optimization method: n n 2 P S H j T ? ? = (5) To determine the SG parameters different optimization functions from Matlab library can be used [4,5,6]: fmincon, fminbnd, fminsearch, fminunc. The nominal time of the machine acceleration Tj is determined by the equation [1]: VII. Comparative Analysis of the SG Parameters Determining Methods At the same time fminunc function provides such methods of smooth unconstrained optimization: ? Steepest Descent method; ? BFGS-method; ? DFP-method (Davidon-Fletcher-Powell). Most simple to use is a function of multivariate unconstrained optimization -fminsearch. In this paper, this optimization function is used. The rest of these the optimization algorithms require a more detailed analysis of their work, setting limits and the initial data. The values of errors of the SG parameters determination by the methods, provided [1] and the optimization method, are shown in Table 1. On excitation winding the test signals are served and, at the same time, the stator phase currents and voltages are fixed. Knowing the input signal and the output signal using the optimization algorithms the SG parameters can be determined. # One of the important directions of the optimization algorithms usage is the synchronous # Global The development of the theme of the SG parameters determination by using the measurements obtained from the PMU can serve as the basis for forming the system of monitoring the synchronous machine parameters condition. It is possible to form the function, where the AEC parameters act as unknown, by using methods of spectrum analysis for PMU measurements processing and analyzing frequency methods of the automatic control theory [16]. Optimization algorithms allow processing this function and finding unknown coefficients. This article deals with the technique of determination of the AEC parameters by using synchronized pharos measurements. # IX. BRIEF DESCRIPTION OF THE INVESTIGATIONAL AEC Development of the technique is based on controller operation analysis of AEC with two stabilization systems (ARV-2SS) of the Russian production. Fig. 4 shows mathematical model of microprocessor controller ARV-2S algorithm, all the model elements are described by the transfer functions of the continuous operator s. As input parameters of stabilization channels in ARV-2SS can be used frequency deviation of voltage of synchronous generator stator and its first derivative, first derivative voltage of generator stator and first derivative of the rotor current [17]. The technique is based on determination of main AEC parameters by optimization algorithms using PMU measurements. # VIII. SG AND AEC PARAMETERS MONITORING BY THE OPTIMIZATION ALGORITHMS USAGE By using the simulation in the hardware and software system RTDS (Real-Time Digital Simulator) [18], PMU measurements have been obtained for different modes of AEC operation. The PMU measurements have Year 2017 F Monitoring of the Synchronous Machines Parameters and their AEC using Synchronized Vector Measurements been formed in data sets, which have been further processed in MATLAB using optimization algorithms. These algorithms find minimum of the objective function, which represents the difference between transfer function obtained from PMU measurements and transfer function calculated with the unknown coefficients of the AEC channels. During technique development, several optimization algorithms have been considered, they are described below. All algorithms use the same objective function. # a) Fminunc algorithm The function Fminunc(fun, x0, options) from Optimization Toolbox MATLAB implements unconstrained nonlinear optimization methods [19]. Each experiment is calculated with ten starting points (x0), chosen at random within the possible coefficients range. Fig. 5 shows the process of calculating coefficients K0u and Kint. Different initial values (starting points) are marked with colours. The convergence is observed for all initial values after about 150 iterations. Main coefficients of AEC channels are listed below. These coefficients are used in this technique as unknown. X. The function Fmincon from Optimization Toolbox MATLAB attempts to find a constrained minimum of a scalar function of several variables so this function implements constrained nonlinear optimization or nonlinear programming [20]. The algorithm allows to restrict the range of variation of the unknown coefficients. The results obtained by using this algorithm fully match with the results of calculation by Fminunc algorithm. However, it must be noted that in some cases this algorithm calculates faster. # DESCRIPTION OF tHE DEVELOPED TECHNIQUE # c) Genetic algorithm The possibility of using genetic algorithm (GA) to calculate the AEC parameters have been considered. The genetic algorithm is a method for solving optimization problems based on biological principles of natural selection and evolution. In this paper, we have applied a standard GA without modification, which is presented in MATLAB. Since the GA uses the principles of calculation different from that of the Fminunc algorithm there are differences in the presentation of the results (mismatch of scale on the horizontal axis in Fig. 5 and Fig. 6). However, the calculation results of both algorithms coincide. The calculation time for any presented algorithms is considerably increased when increasing the number of unknown parameters in the objective function. # d) Comparison of optimization algorithms During the comparative analysis of considered algorithms, the calculation time and maximum errors of calculation have been estimated. As a result, despite the good performance and sufficient accuracy, all things being equal, the GA requires much more time to calculate than the other algorithms, especially when increasing the number of unknown parameters. Therefore, further calculations have been performed using Fminunc algorithm. ARV channels work, and all the channels of system stabilizer are turned on Fig. 8 shows amplitude spectrums of oscillations of the rotor voltage (Vf) for each modes of AEC operation. The number under the picture corresponds to the number of operating mode (Table II) for which it has been plotted. The amplitude spectrums of oscillations of the rotor voltage, based on PMU measurements, are plotted in figure by the solid green curve. The amplitude spectrums of oscillations of the rotor voltage, obtained through calculated coefficients, are plotted in figure by the dotted red curve. In all cases, the curves coincide with sufficient accuracy. The presented technique uses two methods to form the objective function for determining of the AEC parameters by using obtained measurements. The first method uses transformation from spectrum of oscillation of the rotor voltage Vf to spectrum of oscillation relevant input signals of stabilization channels (Ug, fg, If) in the closed automatic control system (ACS), while it is necessary to know a transfer functions of synchronous generator SM through relevant variables. When these data on a synchronous generator absence, the objective function forms by transformation from the spectrum of oscillation input signals of stabilization channels (Ug, fg, If) to spectrum of oscillation of the rotor voltage Vf in opened ACS. This is the second method to form the objective function. The results of calculation of the AVR channels parameters in modes, listed in Table 2, are shown in Fig. 9. The results of calculation of the coefficients of internal and external stabilization channels (Fig. 10) are shows only for modes in which these coefficients are non-zero. # XI. THE RESULTS OF USING THE TECHNIQUE Each figure shows the calculation results for both methods. The errors of the coefficients calculation are mostly less than 1%. In some cases, the error can reach 3%. The second method allows to perform calculations quickly. However, when calculating the coefficients Kint, K1if, K1u and K0f, there is an increase of the error unlike first method. Also in the study the comparative analysis of applying different optimization algorithms has been performed. As a result, despite the good performance and sufficient accuracy the GA requires much more time to calculate than the other algorithms, especially when increasing the number of unknown parameters. # Global The developed technique determines AEC parameters with sufficient accuracy by using PMU measurements, including when in the circuit there is a noise of different power level. Thus, the technique allows to extend the field of PMU application, and with further development, will allow to monitor AEC and calculate its parameters in real time. # Global ![Monitoring of the Synchronous Machines Parameters and their AEC using Synchronized Vector Measurements T.G. Klimova ? , A.I. Kovalenko ? & O.O. Nikolaeva ?](image-2.png "") 1![Figure 1: Saturated open-circuit characteristic (a), unsaturated, (c), the characteristic of the three-phase short-circuit (b) built for the SG model](image-3.png "Figure 1 :") 2![Figure 2: The amount of the Subtransient and transition phase current components (a), a transition component (b), subtransient component (c) of the short-circuit current](image-4.png "Figure 2 :") 72017![Journal of Researches in Engineering ( ) Volume XVII Issue VI Version I YearMonitoring of the Synchronous Machines Parameters and their AEC using Synchronized Vector Measurements V.](image-5.png "7 2017 F") 82017![nominal speed; P n -nominal active power; S n -nominal apparent power; P mech -mechanical losses at the nominal speed; P steel -iron losses at a voltage of experience and the nominal speed. Storage-energy constant H can be determined by the equation[1]:© 2017 Global Journals Inc. (US)Global Journal of Researches in Engineering ( ) Volume XVII Issue VI Version I YearMonitoring of the Synchronous Machines Parameters and their AEC using Synchronized Vector Measurements](image-6.png "8 2017 F") 92017![Journal of Researches in Engineering ( ) Volume XVII Issue VI Version I YearMonitoring of the Synchronous Machines Parameters and their AEC using Synchronized Vector Measurements machine parameters determination in real-time environment or at specified sampling rate. Differential equation system of the Park-Gorev can be used to solve this problem. All the phase stator currents and voltages and rotor current and voltage are recorded by means of phasor measurement units (PMU) and are known quantities.](image-7.png "9 2017 F") 4![Figure 4: ?athematical model AEC with two stabilization systems AFV-2SS The parameters of the automatic voltage regulator (AVR) channel: ? K0u, pu -proportion gain on voltage deviation of generator stator; ? Kint, sec -integral gain in circuit regulating the stator voltage. ? The parameters of the external stabilization: ? K0f, pu/Hz -gain on frequency deviation; ? K1f, pu/Hz/sec -gain on the first derivative frequency. ? The parameters of the internal stabilization:K1u, pu/sec -gain on derivative of the stator voltage; ? K1if, pu/sec -gain on the first derivative of the rotor current.](image-8.png "Figure 4 :") 5![Figure 5: Calculation of coefficients K0u (a) and Kint (b) by Fminunc algorithm b) Fmincon algorithmThe function Fmincon from Optimization Toolbox MATLAB attempts to find a constrained minimum of a scalar function of several variables so this function implements constrained nonlinear optimization or nonlinear programming[20]. The algorithm allows to restrict the range of variation of the unknown coefficients. The results obtained by using this algorithm fully match with the results of calculation by Fminunc algorithm. However, it must be noted that in some cases this algorithm calculates faster.](image-9.png "Figure 5 :") 62017![shows process of calculating coefficients of AVR channel. The horizontal axis denotes the number of individuals, i.e. the number of possible solutions of algorithm.Global Journal of Researches in Engineering ( ) Volume XVII Issue VI Version I 11 YearMonitoring of the Synchronous Machines Parameters and their AEC using Synchronized Vector Measurements](image-10.png "Fig. 6 2017 F") 6![Figure 6: Calculation of coefficients K0u (a) and Kint (b) by genetic algorithm](image-11.png "Figure 6 :") 7![Figure 7: Column diagram of calculation time (a) and maximum errors of calculations (b) for considered algorithms when all stabilization channels are turned off](image-12.png "Figure 7 :") 8![Figure 8: Amplitude spectrums of oscillations of the rotor voltage for different modes of AEC operation](image-13.png "Figure 8 :") 2017![Journal of Researches in Engineering ( ) Volume XVII Issue VI Version I 13 Year Monitoring of the Synchronous Machines Parameters and their AEC using Synchronized Vector Measurements](image-14.png "2017 F") 9![Figure 9: Combined graphs of the values of the coefficient K0u(a) Kint(b) and its error calculation for different modes of AEC operation](image-15.png "Figure 9 :") 1SG parameterThe value set in the SG model (data sheet), p.u.The value obtained by the method of GOST (traditional p.u. method),The value p.u. obtained by the optimiza tion method,Error of the method according to GOST, %Error of the optimization method, %Xd1,711,6921,70981,050,012X'd0,1720,1570,171998,720,006X''d0,1190,1090,118998,40,008T'd0,770,7640,769980,780,003T''d0,09620,1080,9610412,270,1Xl0,180,150,17998116,670,011T'd010,9990,9934370,10,66Rs0,00280,0018970,00280132,310,08Ta0,1320,1610,11921,979,85X 20,1450,11990,130917,39,72H2,1362,02592,1360355,160,0016Tj5,34065,064755,3408755,160,0052 2Mode numberMode description1Only ARV channels work, and the channels of system stabilizer are turned off2ARV and external stabilization channels work3ARV and internal stabilization channels work4 Year 201715F © 2017 Global Journals Inc. (US) Also the presented technique has been checked, when in the circuit there is a noise of different power level. For these experiments, we have used ower system model in RTDS. As previously, the calculations have been carried out in two methods using the Fminunc algorithm with six unknown parameters in objective function. The technique is fully confirmed. ## XII. CONCLUSION In this paper several synchronous generator parameters were determined; the other SG parameters can also be determined by using the optimization algorithms. These results confirm the effectiveness of the proposed alternative methodology for the SM model parameters determining. The error in the parameters determining of this method is much lower than with classical methods and almost absent. The error in determining Ta and X2 parameters involves the usage of the approximate calculation equations. On further research this error can be reduced. Method using the optimization algorithms can also be used in combination with traditional methods. The developed method can also be used to verify the existing SM model, taking into account the fact that the needed oscillograms and characteristics of simulated machine exist. The proposed method is universal and can be used in the different PCs and PACs. Since the parameters are selected by iterations directly in the simulation environment, the method automatically takes into account the peculiarities of simulation SG in any given PCs or PACs. Developing the direction of application of optimization algorithms to determine the SM parameters can achieve even more accurate forming the SM model, which will increase the quality of carried by them * Electric three-phase synchronous machine. Test methods] -Instead of GOST 10169-68 vved. 1978-01-01 Moscow: Izdvo standartov Mashiny jelektricheskie trehfaznye sinhronnye. Metody ispytanij. 1984 85 p * Metody raschetov rezhimov raboty electrooborudovaniya electricheskih stancij i podstancij [The calculation methods of operating modes of electric power stations and substations YPKuznecov * Electricheskie mashiny. [Electric machines] The textbook for university students 1978 Energy Publ 3 edition., revised -L.. 832 p * Ispol'zovanie paketa Matlab v nauchnoi i uchebnoi rabote NYZolotyh Use of MATLAB in scientific and educational work * Informatsionnie tehnologii i komp'youternay matematika" [Educational materials for the training program 2006 Nizhny Novgorod Information technology and computer mathematics * AFDaschenko VHKirillov LVKolomiec VFOrobei Matlab V Inzhenernyh I Nauchnyh Raschetah Use of MATLAB in the engineering and scientific calculations * AstroprintOdessa Publ 2003 * Modelirovanie elektrotehnicheskih ustrojstv v MATLAB, SimPowerSystems i Simulink [Simulation of electrotechnical devices in Matlab IVChernyh SimPowerSystems and Simulink * SPB: Piter 288 2008 DMK Press * Definitions and update of synchronous machine model parameters using optimization algorithms // Sbornik materialov ???VIII sessii Vserosiiskogo nauchnogo seminara po tematike «Diagnostika energooborudovaniya AIKovalenko TGKlimova 2016 Novocherkassk Diagnostic of power equipment * Application of optimization algorithms for definitions and update of synchronous machine model parameters // Materialy VII Mezhdunarodnoi molodezhnoi nauchno-tehnichskoi konferencii «Electroenergetica glazami molodezhi -2016 AIKovalenko TGKlimova 2016 Kazan Power engineering through the eyes of the youth -2016 * Ob'em I Normy Ispytanij Jelektrooborudovanija RD 34.45-51.300-97 Shestoe izdanie. Utv. RAO «EJeS Rossii» 08.05.1997. M., Izdatel'stvo NC JeNAS 2004 * Vvedenie v MATLAB LAMironovskij KPetrova Ju SPbGUAP. SPb 2005 Introduction to MATLAB] Uchebnoe posobie * Jelektricheskie mashiny: Uchebnik dlja vuzov IPKopylov Jenergoatomizdat 1986 360 s.: il * Trebovaniya k sistemam vozbuzhdeniya I avtomaticheskim regulyatoram silnogo deistviya [Requirements for excitation systems and automatic regulators of excitation of strong action C37Ieee Std 13. STO 59012820.29.160.20.001-2012 118 1-2011 -IEEE Standard for Synchrophasor Measurements for Power Systems * Electromagnitnie perehodnie processi v electricheskih sistemah [Electromagnetic transients in electrical systems SA ; MUl'yanov Energiya 1970 520 Uchebnik dlya energeticheskih vuzov i fakultetov * Perehodnie Electromehanicheskie processi ? electricheskih sistemah VAVenikov Ucheb. Dlya electroenerget. Spec. Vuzov. -4 izd., pererab. i dop. -?.: Vissh. Shk. 1985 Transient electromechanical processes in electrical systems. 536 p * Fundamentals of the theory of automatic control. Frequency methods of analysis and synthesis systems AENikulin ??? 59012820.29.160.20.001-2012, approval date 03.04.2012 2004. 2012 17 Saint-Petersburg Requirements To Excitation Systems And Automatic Pid Regulators Of Synchronous Generators * Application of software and hardware system RTDS for the analysis of automatic excitation controller: obtaining and verification of the models of microprocessor AECs JLArtsyshevsky TGKlimova AVZhukov EISatsuk AIRasscheplyev 2014 1 Energetic, Moscow * Using Matlab package in the scientific and educational work. Educational materials for training program NYuZolotyh 2006 Nizhnij Novgorod * Matlab in the engineering and scientific calculations OFDashchenko VHKirillov LVKolomiec VFOrobej 2003 Astroprint, Odessa