processor SYNTHESIS AND TECHNICAL REALIZATION OF CONTROL SYSTEMS WITH DISCRETE FRACTIONAL INTEGRAL-DIFFERENTIATING

Control systems with a fractional order which provide better dynamic and static indicators for many technical objects in comparison with systems with integer order of astaticism were studied. Based on the analysis of frequency characteristics, transient processes and a modified criterion for quality assessment, optimal relationships between parameters of the desired transfer function were obtained. Normalized transition functions of closed systems with the order of astaticism from 1 to 2 were presented with overregulation less than 2...5 % on the basis of which parameters can be chosen and the controller structure determined. The process of stabilizing cutting power was analyzed for a milling machine as an example of the systems with nonlinear parametric and structural dependences in control and perturbation channels. It was shown that fractional integral-differentiating controllers make it possible to provide the order of astaticism from 1.3 to 1.7 and permissible level of overregulation in a wide range of external perturbing influences. A method for approximate calculation of fractional integrals based on the approximation of the highest coefficients of expansion in a series of geometric progressions was developed. It provides reduction of the memory capacity required to store the coefficient arrays and the history of the input signal and requires significantly less CPU time to calculate the controller signal. For example, for controllers based on the Intel® Quark ™ SoC X1000 or FPGA Altera Cyclone V, the quantization period is 6...15 μs and several milliseconds for Atmega328. This makes it possible to implement fractional integral-differentiating controllers based on widely used modern processors and apply fractional-integral calculus methods for synthesis of high-speed automatic control systems. The proposed methods can be used in the control of the objects both with fractional and integer orders of differential equations.


Introduction
The fractional-integral calculus extends the theory of differential equations to the domain of non-integral order of derivatives eliminating discontinuities of this parameter. Solutions of linear equations of this kind represent a wide class of various functions that differ from the combination of exponential and harmonic functions. These functions are the infinite oscillating series. Recently, from 1993 to 2007, when developing the theory of Mittag-Leffler functions, solutions to some fractional differential equations were obtained in general form. These functions are named after their developers: Robotnov-Hartley, Erdei, Miller-Ross. Unification of these functions in 2007 has led to the development of generalized R-and G-functions. Despite the complexity of these functions, fields of application of the fractional-integral calculus in modern technical systems are gradually expanding. Creation of the chaos theory [1] has combined mathematical description of fractal spatial and surface structures and interaction of mobile particles, which is rather effective in describing the processes occurring in porous structures of filters, catalysts, accumulator electrodes and supercapac-itors [2]. Application of the methods of fractional-integral calculus to the wide list of physical processes given in [3] enables obtaining the best results.
One of the features of solution of fractional-differential equations is their close connection with power functions. Accordingly, the control objects in a number of technological processes with nonlinear dependences approximated by these functions are better described by fractional-differential equations [3]. Currently, methods for identifying such objects by means of differential equations with an arbitrary fractional order are already well tested. Accordingly, to control this kind of processes, it is advisable to use PI α D β -controllers with fractional-integral and fractional-differentiating components, which ensure setting of the closed loop to obtain the required dynamic and static parameters.
The theory of synthesis of fractional PID controllers has been well developed [4]. However, the problem consists in technical implementation of such controllers since the use of discrete Grünwald-Letnikov or Riemann-Liouville forms which are infinite rows theoretically involves allocation of infinite memory capacities and requires a large number of arithmetic operations to be performed during the processor quantization period. This limits the scope of application of the PI α D β controllers controlling slow process. The issues of choice of coefficients of the PI α D β -controllers when working with non-linear objects with variable parameters remain open as well.
Therefore, the problems of synthesis of optimized systems with a fractional order of astaticism for objects with variable parameters and development of methods for technical implementation of fast-acting fractional integral-differentiating controllers remain topical.

Literature review and problem statement
Effectiveness of application of fractional-integral calculus is no longer open to doubt in many areas of science and technology. In the modern theory of heat and mass exchange, such a mathematical apparatus enables more accurate solutions for the physics of fractal media in comparison with the methods of integer integration [5][6][7]. For example, it was shown in [6] that dependence of temperature on time in a thermal system operating in the transient mode can be approximated and modeled by means of a differential equation of fractional order. For some particular cases, analytical solutions based on the Mittag-Lefler functions were obtained and the use of PID controllers of non-integer order for process control was proposed [8]. However, in a number of cases, the structure of controllers should be more complex and standard methods for tuning PI α D β -controllers may be ineffective.
Application of the fractional calculus in automatic control can be conditionally divided into two fields. One of them includes the methods of mathematical and computer simulation of the fractional order systems in which properties of fractional dynamics manifest themselves. Discrete calculation methods with a constant step are used there for calculating fractional integral and differential signal components.
The Grunwald-Letnikov (GL) form is a consequence of generalization to an arbitrary order of the Cauchy formula of repeated differentiation and determines the rules for computing fractional derivatives: The complexity of application (1) is caused by the necessity of calculation of the gamma function of large arguments and the rapid loss of calculation accuracy because of the disparate values of the coefficients.
In its turn, the discrete Riemann-Liouville (RL) form is a generalization to an arbitrary order of the Cauchy formula of iterated integration and is used to calculate the fractional integral: The error of this formula is caused by the presence of the integrand denominator which becomes equal to 0 at t=t i . If the last summand is neglected (which is usual for the method of finite increments), then the error is very noticeable, especially at the beginning of the transient process. There is a method for calculating the correcting function [4] calculated from formula however, when calculating the correction value, it is necessary to calculate the gamma function and the factorial with arguments equal to the number of the calculated points which (as already emphasized) leads to certain problems. Besides, calculation of correction is algorithmically equivalent to calculation of one more fractional integral. The second field includes the methods using the apparatus of fractional calculus for synthesis of the systems controlling the objects of both integer and fraction orders, in particular, synthesis of controllers of fraction order [10]. Because of their transcendence, the resulting transfer functions significantly complicate solution of the classical problems in the theory of automatic control: estimation of stability, observability, robustness [11]. However, after the synthesis of controllers, problems of their implementation arise because of complexity of calculating infinite rows in microcontrollers. Therefore, some approximate analog models of such controllers are proposed based on serial RC-circuits [12] which in a case of sufficiently large number of links become similar to fractals and the number of components is limited in the controllers. Also, the fractional-integral components are replaced by high-order filters which are sufficiently similar in their properties in the mid-frequency range [13].
Metal processing is one of the technology fields where the PI α D β -controllers can provide an improvement in quality. The metal-cutting machines operate in conditions of action of parametric perturbations caused by variation of the depth and width of cutting, the workpiece material hardness and the cutting tool condition. Stabilization of the cutting power in machine tools due to the feed speed control helps to rise machining productivity and precision and extend the cutting tool service life. However, in combination with yielding of the "machine-part" system elements, perturbations cause changes in the equivalent parameters of the control object which results in a deteriorated quality of the transient processes. Therefore, the use of classical control systems with anticipatory correction, flexible feedback in the regulated coordinate does not ensure achievement of a high quality of transient processes and smaller static errors in the steady state at a wide range of changes in processing conditions. For metal rolling mills, surface roughness of the worked billets is characterized in [14] as a fractal. This has made it possible to develop a high-precision system for assessing the product quality. It can be assumed that good results can be obtained for other metal processing methods as well.
On the other hand, development of the microprocessor technology has led in recent years to the mass spread of universal single-chip microprocessors manufactured using the 14-28 nm technology processes. Such processors are already used not only in personal computers and smartphones but also in Arduino, Raspberry Pi and Altera controllers.
Thus, some issues remain open concerning determination of parameters and structure of integral-differentiating controllers of a fractional order. The issues of the possibility of applying the theory of fractional integral calculus to the processes of metal processing by cutting are still open. There are mathematical and technical difficulties in implementing the PI α D β -controllers based on inexpensive microcontrollers. Solution of these problems opens the way for building systems for control of various technological objects on the basis of complex mathematical models and controllers without simplifications which was impossible previously.

The aim and objectives of the study
The study objective is the synthesis of fractional integral-differentiating controllers that provide the required order of astaticism at a small overregulation for the systems with increased requirements to dynamic and static indicators as well as development of a method for implementing high-speed discrete fractional PID controllers based on modern microprocessors.
To achieve this objective, the following tasks should be solved: -to develop a criterion for synthesis of controllers to ensure optimal dynamic and static indicators; -to check operability of regulators for an object with parameters varying in a certain range; -to develop a method for implementing discrete fractional-integral high-speed controllers corresponding to the capabilities of modern single-chip microprocessors.

Synthesis of closed systems with optimal dynamic and static indicators
Let us consider a closed system in which we shall use not a PI α D β -controller but provide setting to a given optimum. High accuracy of stabilization of coordinates in static and dynamic modes can be provided by astatic systems with their open loop containing an integrating link of the first or higher order corresponding to the required order of astaticism.
Setting to the modular optimum with a calculated overregulation of 4.35 % is widely used for the first-order astatic systems and setting to a symmetrical optimum at which overregulation increases by a factor of 10 is used for the second-order astatic systems. The reason for this is reduction in phase and amplitude stability margin.
To unify the subsequent calculations, move to relative units in the time domain. This will allow us to consider from a single point of view both relatively slow systems (e. g. climate systems) and high-speed systems (such as metal cutting machines, welding machines). Assume where t, t real are relative time and absolute time, T ν is uncompensated small time constant of the control object. By selecting this value, one can set the required system speed based on the requirements to the process. In particular, a number of additional restrictions related to equipment wear and tear, financial and energy costs may be taken into account [18].
Then, the Laplace operator and frequency will also be changed in the transfer functions and in calculation of frequency characteristics: Synthesize a system with a fractional order of astaticism and the open loop transfer function described by the following expression: where a, b are setting parameters and the 1 1 p + component corresponds to the uncompensated part of the control object. Such a system does not have a position error. A constant velocity error arises when a "concave" signal of the t 1+µ form is applied to the input.
The transfer function (6) corresponds to the logarithmic amplitude-and phase-frequency characteristics (LAFC and LPFC) calculated from formulas: Comparison of LAFC and LPFC of an open loop with a fractional order of astaticism and the systems with tuning to a modular or symmetric optimum at identical cutoff frequencies shows that the phase stability margin increases by ( ) with a decrease in µ. This makes it possible to improve dynamic performance of the system by increasing the cutoff frequency while maintaining or even increasing the stability margins. The transition function (6) of the system which is a convolution of the Robotnov-Hartley function and damped harmonic oscillations [4] is characterized by several extrema. For each µ, one can choose a certain a to b ratio that satisfies the set quality criterion of the transient process.
The classical criterion of the mean-square total error can be formally represented as a sum of the products of multiplication of errors by their "weight" equal to the error itself. As a result, the largest "contribution" is made by the first points with the maximum deviation and the process "pressed" to the ordinate axis appears to be optimal. The overregulation can be very large but short-lived.
Let us formulate a criterion for optimizing the transient process for the systems requiring provision of high system speed (comparable to the chosen value of T ν ) but with a strict limitation of width of the admissible "corridor" of deviations from the set value. To this end, integrate the error before the first entrance to the ( ) 1 i Y − < δ "corridor" (the error weight is 1) and then increase the F functional by the amount of the cubed error relative to δ. Accordingly, the weight appears to be less than 1 in the "corridor" and more than 1 outside it and rises according to a quadratic parabola: where t 1 is the time of the first entry into the "corridor". This formulation of the estimation criterion makes it possible to avoid compression of the optimal transition function to 0 in time and short-term surges and depressions outside the set corridor [15]. Fig. 1, 2 show the normalized transition functions obtained based on the set criterion. Optimum characteristic can be chosen from them on the basis of the desired system speed. The corresponding combinations of the µ, a, b values can be taken from Table 1.
The values of the margin of Δφ phase stability and Ω 0 cutoff frequency for δ=0.025 and δ=0.05 are also given in Table 1. It can be seen that an increase in δ leads to a certain decrease in the phase stability margin with a simultaneous increase in the cutoff frequency, that is, to an increase in the system speed. Also, for all combinations of the parameters obtained, the phase stability margins are even greater than when tuned to the modular optimum (65.5°) and overregulation is smaller. Table 1 Optimal ratios of the setting parameters

Analysis of operation of the cutting power stabilization system of the milling machine with fractional-integral controllers
Let us check how the system behavior changes with the change of the control object parameters. Consider a system for stabilizing the milling machine cutting power as such an object. The feed rate is the control coordinate in this machine and a sudden change in the cutting depth is the perturbation action.
The static value of the cutting power is determined from the empirical relationship [16]: where K is the equivalent coefficient depending on the types and parameters of the tool and the part, t p is the cutting depth, s is the feed rate,  Fig. 3 shows two of the oscillograms obtained in the study of a vertical milling machine, model 6B75, for milling gray cast iron with an end milling cutter made of R6M5 grade high-speed steel. The active power of the asynchronous motor of the main motion when the cutter cuts into the workpiece was determined by the K-50 measuring instrument of 0.5 accuracy class with further signal filtering by a first-order inertial link with a time constant T F =0.2 s. The "machine-detail" system was identified by classical methods for the transfer function of the form: where A is the coefficient determining condition of the cutting tool and the material of the workpiece, z is the number of teeth of the milling cutter and n is the rotational speed of the milling cutter [rpm]. The synthesized control system with a fuzzy controller has made it possible to reduce cutting power to 15...25 % in one second [17]. Taking into account formula (10), it is obvious that the mathematical description of the control object includes power dependences which, under the Laplace transform, leads to a fractional-differential equation in the operator form. The features of the objects described by such equations are clearly visible in Fig. 3: the rapid onset and the prolonged termination of the transient process.
Then identification of the object by the transient processes from Fig. 3 leads to the results presented in Table 2. Table 2 Results of control object identification  With the relative rms error of F to 1.4...2.1 % determined by noise and see-saw pulsations, the cutting process is described by fractional-differential equations. The order of the equations is comparable with the order of the power dependence (10) for the workpiece and cutter mate- The smoothed lines in Fig. 3 convince that the calculated transient processes are close to the experimental ones.
To ensure the selected settings, the controller connected in series with the control object must have a transfer function determined from the relation Then, In a particular case, when 1, O µ = that is for the objects with an integer order of differential equations, the following is obtained: p structure according to formula (16) for the worst case of object parameters: the maximum equivalent gain and the largest coefficients of the denominator of the transfer function corresponding to the nominal power. Fig. 4 shows response of the cutting power stabilization system to the stepped changes in the cutting depth resulting in cutting power jumps in the open system from 100 % to 140 %, 180 %, 220 %, 160 % and again to 100 %. The model takes into account the power dependences of the cutting power both on the feed and the depth of cutting, the noise and the sensor lag. The controller settings were made for µ=0.5. Similar results were obtained for the order of astaticism 0.3 0.7.

≤ µ ≤
At lower values of µ, the system is sensitive to noise. To make its operation stable, a quantization period ( ) 0.005...0.02 , t T ν ∆ < is required. At higher values of µ, the system speed is much lower than the variants shown in Fig. 4.

Fig. 4. Graphs of transient processes in the system of power stabilization at jumps in the depth of cutting
Based on the obtained graphs, it can be concluded that the proposed method of setting is advisable. Obviously, the synthesized controller provides the desired quality of transient processes: the cutting power in any of the modes does not deviate more than 2 % from the set value.

The method of technical realization of fractional integral calculation for microprocessor-based control systems
To realize calculation of the fractional-integral and differential components of the controller signal in the microprocessor systems with quantization period Δt, several calculation methods can be used. The Grunwald-Letnikov (1) and Riemann-Liouville (2) forms lead to close results. However, considering, that the systems' astaticism provides integral components of controllers and the differentiating components appear only when the inertial components of the object of the 2nd and higher order are compensated, it is more often necessary to perform an integration operation to find solutions. Therefore, we shall prefer the form (2) and use the following relation to calculate the fractional derivative: To compensate for the calculation error by the Riemann-Liouville form, introduce a corrective bias into the denominator. This will allow us to obtain an exact solution for some particular cases and more accurate solutions for arbitrary-shaped signals.
If a single jump is applied to the input of a fractional-integrating link, solution of the equation ( ) ( ) ( ) I 1 y t t µ = is described by expression: Then, write Γ + µ − ∆ + ∑ and the following formula for calculating the fractional integral is obtained after a series of transformations: where j k µ are constant coefficients. The obtained solution (21) was tested using the relation by integrating a pseudo-random signal (of "white" noise).
In the proposed method (21), because of the shift in the denominator, a regular error inherent in the basic form was eliminated without applying corrective calculations. However, the main problem of calculating the fractional integral using either of formulas (1), (2) or (21) consists in a necessity of calculation of sums of the pairwise products of the coefficients by all previous values of the input signal in each quantization period. In this case, the processor memory should store both the array of coefficients and the entire "history" of the input signal. Neglect of the highest terms of the series leads to an error in calculation of ε. It is possible to determine the minimum number of points from (21) to provide the set error value: For example, for µ=0.5 and ε=0.05 at Δt=1, 315 previous points should be remembered and 1,964 points to be remembered to ensure ε=0.02. Storage of this amount of data exceeds capabilities of many single-chip processors. Now, calculate the coefficients for various µ from 0.1 to 0.9 and plot their variation depending on the ordinal number in the logarithmic coordinate axes. At the same time, construct trend lines for these dependences in the form of power functions (Fig. 5).

Fig. 5. Dependence of coefficients in the fractional integration formulas on the ordinal number
A remarkable feature of these dependences is obvious. Starting from 50...100 numbers, the differences between the forms (2), (21) practically disappear and the dependences of the coefficients with an error in the 5 th or the 6 th decimal place are approximated by the power function, that is, they become close to a geometric progression with a coefficient q extremely close to but still less than 1. This means that calculation of the sum of products of higher terms can be replaced by multiplying the previous value of this sum by the coefficient of geometric progression. Then, by setting the maximum number of the stored values of the input signal n dim , calculation of the fractional integral can be performed in two blocks. The initial n dim points are calculated by (21), and then all leading summands are gradually taken into account in the sum of the geometric progression S. Only the last n dim points are processed by (21) and form the dim I .
where q is calculated as the mean geometric between Efficiency of the algorithm was verified by calculating the fractional integral of the sinusoidal signal ( ) ( ) sin x t t = for various µ. Fig. 6 shows a quarter of the period of the results of integration by (24) from the 900 th to the 1,000 th calculation points at a set n dim =1000 (the time axis is marked in degrees). The phase shifts are exactly equal to , 2 π µ that corresponds to the definition of the fractional integral of the harmonic function. The error in the calculation results in comparison with the calculation by (21) is from 0.11 % for µ=0.9 to 1.7 % for µ=0.1. The speed of the calculation method was tested on Arduino Intel Galileo Gen 2 microcontrollers based on the 8086 processor with clock speed of 400 MHz, FPGA Altera Cyclone V with clock speed of 50 MHz in the DE1-SoC controller as well as Arduino Nano with Atmega328 processor with clock speed of 8 MHz. ).
double was used. The exact data and the data calculated by the controller coincide at n max which can serve as a criterion for n max choice. For comparison, the graph of y_old change illustrates the result of calculating the fractional integral when memory capacities are exhausted: the fractional-integrating link becomes proportional and a static error appears in the astatic system. The calculation time for the fractional integral was on average 15 µs, periodic communication procedures took up to 15 µs as well. Calculation time for the fractional integral was 8 µs at n dim =64. FPGA Altera used integer mathematics with data storage in 64-bit registers. Due to organization of parallel computations available in FPGA, the calculation time for the fractional integral was 6 µs. Even for the simplest Atmega328 based controller, the calculation time was 2.2 ms at n dim =128 and 1.1 ms at n dim =64.
Despite the calculation error from percent fractions to several percent, the most important feature of the developed method consists in reduction of memory capacity requirement and the number of operations necessary to implement fractional-integral computations. Fig. 8 shows the graphs of transients in the Arduino Galileo control of the second-order inertial link with a unit gain and sending to the meander input 1500ó2500 The controller corresponds to transfer function (18) at a 1 =0, a 2 =0 and the quantization period of 50 µs. Comparison of the control signal Ur and the output signal Yo confirms operation of the fractional-integral component: the controller signal does not coincide with the feedback signal in quasi-steady modes.
Thus, by modifying the algorithm of calculation of the fractional integral, a possibility of forming control systems with discrete fractional-integral controllers with a quantization period of units to tens of microseconds was achieved. The obtained results allow us to make a conclusion about effectiveness of application of fractional integral-differentiating controllers. Due to the new variable parameter, namely the fractional order of astaticism, it is possible to obtain the combinations of speed, overregulation and speed error in closed systems that are unattainable with the systems of an integer order. Although such controllers are parametric, they can be applied in nonlinear systems, for example, in metal processing machines. It is possible to select such coefficients and the controller order for nonlinear objects which will provide the desired quality indicators in certain ranges of parameter variation.
The complexity of implementing fractional controllers in digital systems is explained by the necessity of calculation of sums of infinite series. Therefore, in the opinion of the authors, the proposed method for approximating the terms of expansion by geometric progression is an essential step in solving this problem and enables creation of high-speed microprocessor-based control systems. This opens up new fields of application of the methods of the fractional-integral calculus. For example, the use of fractional integral-differentiating controllers can become possible not only for theoretical studies and control of slow processes (such as climate systems) but also for a number of high-speed control systems for metal-cutting machines and electric drives of various types.
The issues of rational choice of the quantization period, the controller structure and coefficients for the optimal systems of asynchronous motor and motors of valve control remain open in this area.

Conclusions
1. The control systems with a fractional order of astaticism which provide better dynamic and static indicators for many technical objects in comparison with the systems with an integer order were investigated. Based on analysis of the frequency characteristics, the transient processes and the modified quality assessment criterion (9), optimal relationships between the parameters of the desired transfer function were obtained. Normalized transition functions of closed systems with an order of astaticism from 1 to 2 and overregulation less than 2...5 % were presented. Parameters can be chosen and the controller structure determined on their basis.
2. Analysis of the process of stabilizing cutting power of milling machines as an example of the systems with nonlinear parametric and structural dependences in the control and perturbation channels was performed. It was shown that fractional integral-differentiating controllers make it possible to provide the order of astaticism from 1.3 to 1.7 and the permissible level of overregulation in a wide range of external perturbing influences. This indicates advisability of applying the methods of fractional-integral calculus for synthesis of systems of metal processing by cutting.
3. A method of approximate calculation of fractional integrals based on approximation of the highest coefficients of series expansion of geometric progressions was developed. It enables reduction of the memory capacity required to store the coefficient arrays and the history of the input signal and requires significantly less CPU time to calculate the controller signals. For example, for the controllers based on Intel® Quark ™ SoC X1000 or FPGA Altera Cyclone V, the quantization period is 6...15 µs and units of millisecond for Atmega328. This makes it possible to implement fractional integral-differentiating controllers based on widely used modern processors and apply the fractional-integral calculus methods in synthesis of high-speed automatic control systems. The proposed methods can be used to control the objects both with fractional and integer orders of differential equations.