CC BY
8
-□ □-
Дослгджено замкнет системии з дробовим порядком аста-тизму, як для багатьох техтчних об'ектгв забезпечують кращi динамiчнi i статичт показники поргвняно з системами з щлочисельним порядком. На пiдставi аналiзу частотних характеристик, перехгдних процетв i модифжованого кри-терхю ощнки якостi отримано оптимальт спiввiдношення мiж параметрами бажаног передавальног функщг. Надано нор-моват перехгдт функщг, на пiдставi яких може бути обрано бажаний порядок астатизму системи i визначено структуру i параметри регулятора.
Розглянуто процес стабт^ацп потужностгрЬання фрезерного верстата, як приклад систем з нелтшними параметрич-ними i структурними залежностями по каналах управлтня i збурення. Показано, що дробовi iнтегрально-диференщйтрегу-лятори дозволяють забезпечити порядок астатизму вiд 1.3 до 1.7 та припустиме перерегулювання в широкому дiапазонi зов-тштх впливiв.
Розроблено метод наближеного розрахунку дробових гнте-гралiв, заснований на апроксимащг старших коефiщентiв роз-кладання в ряд геометричною прогреЫею. Це забезпечуе ско-рочення необхгдних обсягiв пам'ятi для зберыання масивiв коефщ^ентъв i кторп вхiдних сигналiв та вимагае значно мен-ших витрат часу процесора. Для процесорiв Intel® Quark™ SoC X1000 та FPGA Altera Cyclone V перюд квантування склав одинищ микросекунд, для Atmega328 - одинищ мЫшекунд. Це дозволяе реалЬувати дробовi iнтегрально-диференцiйнi регу-лятори на базi широко поширених сучасних мiкропроцесорiв i застосовувати методи дробово-ттегрального числення для синтезу швидкодхючих систем автоматичного управлтня. Запропоноват методи синтезу i реалiзащг дискретних дробових ттегрально-диференщйних регуляторiв можуть бути застосоват для об'ектiв управлтня як з дробовим, так i щлочисельним порядком диференщйних рiвнянь
Ключовi слова: астатична система, дробове ттегруван-ня, дробове диференщювання, алгоритм швидкого розрахунку
дробового ттегралу -□ □-
UDC 644.1+004,9:517.9
I DOI: 10.15587/1729-4061.2018.1398921
SYNTHESIS AND TECHNICAL REALIZATION OF CONTROL SYSTEMS WITH DISCRETE FRACTIONAL INTEGRAL-DIFFERENTIATING CONTROLLERS
V. Busher
Doctor of Technical Sciences, Professor* Е-mail: victor.v.bousher@gmail.com A. Aldai ri Postgraduate student* Е-mail: tom_ford1986@hotmail.com *Department of electromechanical systems with computer control Odessa National Polytechnic University Shevchenka avе., 1, Odessa, Ukraine, 65044
1. 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 P/D^-con-trollers 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 Grunwald-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 PIaDP controllers controlling slow process. The issues of choice of coefficients of the PIaDP-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.
2. 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-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 PIaD$-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:
D* f (t ) = lim—Y (-1)'--f (t - j At),
J\> ^o At* ff ' r( j + 1)r( * — j + l/1 7
1 1
D*f (t) I At=const "A* Í kGLjf1-
(1)
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:
I* f (t ) = 4t ff-dt = -^lim Y f (jAt) At, K ) r(*)f (t, —1)1—* r(*)Y(ts — jAt)1-* ,
The error of this formula is caused by the presence of the integrand denominator which becomes equal to 0 at t=t;. 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
e(1) =
At *
r(k+ *) r(k-1+ *)
(3)
At* 1
f (t) I At=const ~ r(*) AimS ^ kRLjf1-j+1.
(2)
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 PIa 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 PIaDP-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.
Then, the Laplace operator and frequency will also be changed in the transfer functions and in calculation of frequency characteristics:
P = TvPreal ^ P*=(Tv)* PL,
n = TvQreal ^Q* = (Tv)*Qli.
(5)
Synthesize a system with a fractional order of astaticism 1+*, * e [0,1] and the open loop transfer function described by the following expression:
3. 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.
4. Synthesis of closed systems with optimal dynamic and static indicators
Let us consider a closed system in which we shall use not a PIa 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
t=
T
(4)
where t, treai are relative time and absolute time, Tv 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].
hoopt(p)= bjP+P -+T, 0<* <1,
ap*p p+1 where a, b are setting parameters and the
1
p +1
(6)
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 ti+l 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:
L = —20lg (aQ*) — 20lg (Q) — —20lgV Q2 +1 + 20lgVJQ+1,
9 = - (1+* )— arctg (Q) + arctg (bQ) =
, ,% (bQ — Q)
= —(1+*)—+ arctg---—f-.
y 6 1 + bQ2
(7)
(8)
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 (1 — * )— 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 Tv) 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 (|Yi — 1 < S) "corridor" (the error weight is 1) and then increase the F functional by the amount of the cubed error relative to 5. 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:
quantization period At< (0.005...0.02)TV, at smaller |a. The system speed becomes relatively low at larger ^ (tx > 3TV).
f=2
(1 - Y )v t < t, ||1 - YI j v t > t
^ min,
(9)
where t\ 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 ji, a, b values can be taken from Table 1.
The values of the margin of A9 phase stability and cutoff frequency for 8=0.025 and 8=0.05 are also given in Table 1. It can be seen that an increase in 8 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
0=0.025 0=0.005
M a b A9, 0 ^0 M a b A9, 0 ^0
0.1 0.080 1.675 82 15.90 0.1 0.085 5.002 82 40.87
0.2 0.267 2.113 77 5.55 0.2 0.189 3.162 75 10.44
0.3 0.574 2.781 73 3.27 0.3 0.483 3.502 71 4.52
0.4 1.146 4.226 70 2.41 0.4 1.021 4.226 69 2.64
0.5 2.420 6.077 69 1.68 0.5 2.284 6.437 67 1.84
0.6 5.023 9.396 67 1.28 0.6 4.742 9.953 67 1.40
0.7 11.97 13.46 71 0.82 0.7 10.67 14.26 68 0.96
0.8 29.26 20.15 74 0.54 0.8 26.08 22.61 69 0.67
0.9 81.62 40.16 72 0.42 0.9 77.06 47.73 68 0.52
h(t)
0.8
0.6
0.4
0.2
1 / 1 / ft / / / / / /
! / ' '/' / !
•7 ' i : / / / / - |i=0.1 -|I=0.3 - |I=0.6 - [1=0.9
lli / •y r-
77'' / -|i=0.7 |x=0.8
w
Fig. 2. Transient characteristics at the optimal ratio of p, a, b for 8=0.05
5. 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]:
P = KsVp t'
(10)
where K is the equivalent coefficient depending on the types and parameters of the tool and the part, tp is the cutting depth, 5 is the feed rate, Xp, Vp are the exponents that vary depending on the materials of the part and tool:
e [0.75,0.95], e [0.65,0.8].
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 Tp =0.2 s. The "machine-detail" system was identified by classical methods for the transfer function of the form:
Fig. 1. Transient characteristics at the ratio of p, a, b optimal by criterion (9) for 8=0.025(9)
The transition functions at ^e [0.3,0.7], are of greatest interest of all obtained results since it is necessary to use the
Ho ( P)--
(T0p + 1)2 (TFp +1)
(11)
which resulted in the dependence of the milling process time constant on the feed rate 5 [mm/min] and the cutting depth tp [mm] as:
k
0
Industry COntfül SUStëmS
T0 = AzU = Azn-- (tf "s-»5 ),
(12)
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.
Fig. 3. Experimental and calculated graphs of the change in the active power of the main motor at tp=4 mm
Let us take TF as an uncompensated time constant and describe the relationship between relative power and feed by transfer functions of orders 1+| and 2+|:
ho (p)=
1
(ai p1+" + a0 p" +1)( p +1)'
rials and a more complicated H20 (p) or simplified H0 (p) model can be used depending on the required accuracy of calculations.
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
hreg (p) =
Then,
hopt ( p)
ho (p) '
(14)
i ( 1)(a1 PP"0 + a0p"0 +1)
hreg (P) =
app
)= (bp + 1)(a2 ppp"0 + a, pp"0 + ao p"0 +l)
hreg (p)= ,, • (15)
app
When *0 = *, the transfer functions of the controllers are converted into a combination PID +I* (PI) of functions:
HRg ( p )=f ^p+j^+^ 1+-V f bj+1,
^ a a ap) p * ^ a ap)
hreeg (P) =
a2b 2
p +—p +
aa
ab ba0 + a1 a0 1 1 ( b 1 ^ — p + —0-1 + — +—\ — +—
: ap J p" I a ap ,
(16)
When *0 ^ * , the transfer functions of the controllers are more complex and it becomes necessary to calculate two fractional integrals of different orders:
Hreg (p) = (bp +1)
1
Oo+Ol 1+_O±
a ap) p" ap
Ho (p)=
O p2+"+ a, p1+" + ao p" +1)( p +1)'
(13)
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
Parameter
Fig. 3, a
HO ( p )
0.7348
1.4820
2.2372
0.0132
HO ( p )
0.7582
1.6427
2.1110
0.5721
0.0110
Fig. 3, b
HO (p)
0.6652
1.0604
2.3182
0.0208
H ( p)
0.6743
1.1907
2.2372
0.6810
0.0140
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-
HReg (p) = (bp +1)
p
a2 p
a
+£l+Ok. )+±_L )
a ap I p" ap J
(17)
In a particular case, when *0 = 1, that is for the objects with an integer order of differential equations, the following is obtained:
Hreg (p)=4((bp+1)(-p -+-)
p" ( a a ap J
HReg (p)=At (bp+1)
a, p2 a, p . a0 . 1
■ + + a ap
(18)
Simulation of transient processes was performed for controllers with a HR£G (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
1
u
n
n
n
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 At < (0.005...0.02)Tv, is required. At higher values of the system speed is much lower than the variants shown in Fig. 4.
Then, write
qi-H f = d p f _
dt At '
(19)
To compensate for the calculation error by the Rie-mann-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 y(t)= I*(1(t)) is described by expression:
y (t )=
1
r(1+ *)
t *.
(20)
At = -
r(*) 1=1 (ti - jAt + Cj )"* r(1 + *)
'iAt )*,
and the following formula for calculating the fractional integral is obtained after a series of transformations:
I*/(t ) = At
j=1
/-( j - 1f) r(1+ *)
/-
j+1 = At* £ jj
j=1
(21)
where k* are constant coefficients.
The obtained solution (21) was tested using the relation
I* (l1-*( / (t )H / (t )dt
(22)
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.
6. 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 At, 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:
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 pair-wise 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 s. It is possible to determine the minimum number of points from (21) to provide the set error value:
N = |1 r(1 +
e r( *)At*
1/*
(23)
For example, for |=0.5 and s=0.05 at At=1, 315 previous points should be remembered and 1,964 points to be remembered to ensure s=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 5th or the 6th 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 ndim, calculation of the fractional integral can be performed in two blocks. The initial ndim 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 ndim points are processed by (21) and form the I" value. The resulting approximate value of the integral is defined as the sum of S and :
I£( f (t)):
I £ f (t ) = I(At £k £) f
;=i
- j+i
V¿e[1, ndlJ,
5 . = f +, (Atk ) + qS
dim ^ i-1
"dim -1
(24)
I£dim f (t ) = I Kk1) f-j+i ViE["dim, j=1
I£ f (t )= S + Pm f (t)
where q is calculated as the mean geometric between kni and k :
q=
(nmax + 1)£-("max)
("dim + 1)£"(«m„1)'
(25)
Accordingly, only the last ndim values of the input signal and the same number of k ...k coefficients should be stored
1 "dim
in the processor memory.
Efficiency of the algorithm was verified by calculating the fractional integral of the sinusoidal signal x (t) = sin (t) for various Fig. 6 shows a quarter of the period of the results of integration by (24) from the 900th to the 1,000th calculation points at a set ndim=1000 (the time axis is marked
n
in degrees). The phase shifts are exactly equal to , 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 1=0.9 to 1.7 % for |=0.1.
The speed of the calculation method was tested on Ar-duino 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.
Fig. 7 demonstrates comparison of the calculated values (y _ p) and the values (y _ ard) obtained via serial port from Arduino Galileo for a fractional-aperiodic link with |=0.5 at ndim=128, nmax=1000. Double-precision floating-point math (double) was used. The exact data and the data calculated by the controller coincide at nmax which can serve as a criterion for nmax 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.
Fig. 6. Fractional integrals of the sinusoidal function
Fig. 7. Calculated and experimental values of the output signal of the fractional-aperiodic link
The calculation time for the fractional integral was on average 15 |is, periodic communication procedures took up to 15 |is as well. Calculation time for the fractional integral was 8 |s at ndim= 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 |is. Even for the simplest Atmega328 based controller, the calculation time was 2.2 ms at ndim=128 and 1.1 ms at ndim=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 150002500 The controller corresponds to transfer function (18) at a1=0, a2=0 and the quantization period of 50 |is. 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.
Yo, Ur
2000
1000 I
0 " —— — —
0 0.1 0.2 0.3 0.4 t, s
Fig. 8. Graphs of transients in Arduino Galileo
7. Discussion of results obtained in the synthesis and technical realization of discrete fractional integral-differentiating controllers
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.
8. Conclusions
1. The control systems with a fractional order of astat-icism 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 |is 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.
References
1. Mandelbrot B. B. Fraktalyi i haos. Mnozhestvo Mandelbrota i drugie chudesa. [Fractals and chaos. Mandelbrot set and other wonders]. Izhevsk: NIC "Regulyarnaya i haoticheskaya dinamika", 2009. 392 p.
2. Freeborn T. J., Maundy B., Elwakil A. S. Fractional-order models of supercapacitors, batteries and fuel cells: a survey // Materials for Renewable and Sustainable Energy. 2015. Vol. 4, Issue 3. doi: https://doi.org/10.1007/s40243-015-0052-y
3. Uchaikin V. Metod drobnyh proizvodnyh [Method of fractional derivative]. Ul'yanovsk: «Artishok», 2008. 512 p.
4. Das S. Functional Fractional Calculus for System Identification and Controls. Springer-Verlag Berlin Heidelberg, 2008. 240 p. doi: https://doi.org/10.1007/978-3-540-72703-3
5. Tarasov V. E. Fractional Liouville and BBGKI equations // Journal of Physics: Conference Series. 2005. Vol. 7. P. 17-33. doi: https://doi.org/10.1088/1742-6596/7/1/002
6. Svobodnokonvektivnyie techeniya, teplo- i massoobmen [Freeconventional flow, heat and mass transfer] / Gebhart B., Dzhaluri-ya Y., Mahadzhan R., Sammakiya B. Moscow: Mir, 1991. 678 p.
7. Aoki Y., Sen M., Paolucci S. Approximation of transient temperatures in complex geometries using fractional derivatives. Technical Note of department of aerospace of Notre Dam, 2005. 21 p.
8. Petras I. The fractional - order controllers: Methods for their synthesis and application. URL: https://arxiv.org/abs/ math/0004064
9. Jumarie G. Path probability of random fractional systems defined by white noises in coarse-grained time. Application of fractional entropy // Fractional Differential Calculus. 2011. Issue 1. P. 45-87. doi: https://doi.org/10.7153/fdc-01-03
10. Fractional dynamics and control / D. Baleanu, J. A. T. Machado, A. C. J. Luo (Eds.). New York: Springer, 2012. 309 p. doi: https:// doi.org/10.1007/978-1-4614-0457-6
11. Xue D., Chen Y. Q., Atherton D. P. Linear Feedback Control: Analysis and Design with MATLAB. SIAM, 2007. 354 p. doi: https://doi.org/10.1137/1.9780898718621
12. Petras I., Podlubny I., O'Leary P. Analogue realization of Fractional Order Controllers. FBERG, Tech. University of Kosice, 2002. 84 p.
13. Marushchak Y., Kopchak B. Synthesis of automatic control systems by using binomial and Butterworth standard fractional order forms // Computational problems of electrical engineering. 2015. Vol. 5, Issue 2. P. 89-94.
14. Krivolapova L. I., Kravcova O. A., Sokolov S. V. Fraktal'naya razmernost' - ocenochnaya mera kachestva poverhnosti metallo-prokata // Doklady TUSURa. 2015. Issue 1. P. 142-147.
15. Busher V., Yarmolovich V. Modeling and Identification of Systems with Fractional Order Integral and Differential // Electro-technical and Computer Systems. 2014. Issue 15 (91). P. 52-56.
16. Spravochnik po obrabotke metallov rezaniem / Abramov F. N., Kovalenko V. V., Lyubimov V. E. et. al. Kyiv: Tekhnika, 1983. 239 p.
17. Vodichev V. A., Gulyu M. V., Muhammed M. A. Primenenie fazzi-regulyatora v elektromekhanicheskoy sisteme avtomatizacii metalloobrabotki // Visnyk NTU «Kharkivskyi politekhnichnyi instytut ». 2005. Issue 45. P. 504-505.
18. Development of the method for determining optimal parameters of the process of displacement of technological objects / Lutsen-ko I., Tytiuk V., Oksanych I., Rozhnenko Z. // Eastern-European Journal of Enterprise Technologies. 2017. Vol. 6, Issue 2 (90). P. 41-48. doi: https://doi.org/10.15587/1729-4061.2017.116788
