PAR P O M I A RY • A U T O M AT Y K A • R O B O T Y K A
1/2019 ISSN 1427-9126 Indeks 339512
Cena 25,00 zł w tym 5% VAT
Technical Sciences Quarterly | - .
!
W numerze:
3 5
Od Redakcji
The Fractional Order PID Control of the Forced Air Heating System
11
17
! "
21
# " $ % $ #
31
& ' ( $ ) *& * ) '
!
" #
$ % ! % &
'( ( %
( (
# '
'( ) * #
' # (
# ( # &+,
+ Informacje dla Autorów – 59 |
63 | 64 | Awans naukowy ! " " # $ % 65 | 66 | &'" ( ) *$ +,'- 67 | .. ) $ )
+,+, 68 |
% ! ! ! / ! 0 70 | $ $ $ 72 | Polecane
$ !
# # ! 0# *$ 73 | # 1 $ ) +,'- 74
Rada Naukowa Rok 23 (2018) Nr 1(231) ISSN 1427-9126, Indeks 339512
Redaktor naczelny
' %2 3
. # ) > )
JM ' %2 -
, W X Y Z( [
& ; < '
' %2 4 5 K # D K ) \ ]
/ 0'
#
!& = '
' %2 # "2 6 % % C W X Z@ [
/ ' (
7 # + W X > ) Z% ) [
!& > ) '
& ? @ = A # )
& !& = ' A ) B
& # C# A ) =
& ; < ' A $
( 1 D
' %2 ( 8 ( H > J# A ) K # # ) ) M K ? ' %2 9 .2 8 % W X ZW [ ' %2 +# &2 8 ' ? ! ( < ^ ) #
Druk
' %2 32 4 - K # Z #= [
< @ E C # # " ( & & & % FGG = &
' %2 3 - # W X Z( [ ' %2 . - ( % # ) ( H ) # ?
Wydawca ( H > J# A ) K # # ) L ) M K & ) OGO GO PQF ?
Kontakt " # = ) # ) " $ & ) OGO GO PQF ? & OO QTP GU PF # V & & & Pomiary Automatyka Robotyka ) ) # ) $ ) # UccT & UQ # ) & #
#
# ] L
# $
#) # ) $ ) = & # ) # ) " $
L $ > <CD;b @ = ( K%*Df ; D"%K;W( ZK;g
OGUF+ Th OP[
L ! L $ # $ ! )
"K %C & # % (? #$ # Q Z & UOPP[&
] ] % # # ) $
# # #$ L # # )
) # ) " $ & ? ] ] Z ][
&
? ) # ) " $ E ) #) cGTB *W%BOGUc ^ M % # ( ? ! =
^H # ] ] # &
" :2 - W X Z@ [ ' %2 3 $ ? < ] > ' %2 , .2 $ K # ; ( Z" [ ' %2 & _ X X ` % W X ZW [ ' %2 # ! ? K ! < ] ? ! ( D =
< ] ? ' %2 3 ; W X Z [ ' %2 ! ; a W X Z? > [ ' %2 ? K ! # = , #
C> A K #
( ) C = > ) @)$b
Z% ) [ ' %2 # ; < ? C = \ ] ' %2 & = ( X ) ( Z( [ ' %2 8 % 4 < K # K ) # ) " $
? ' %2 ( " KC W X $ # Z # [ ' %2 # 42 " _ X ` % W X ZW [
) # ) " $ K((% UPOT cUOF "& Oh % UBOGUc
3
Od Redakcji
5
The Fractional Order PID Control of the Forced Air Heating System ( K* # ) = # # ) =
11
dynamicznych # *i i ( D ) * ) ( )
17
! "
! " # $ % ) ^ ) # #! ) $
$ M
21
# " $ % $ # &' ' ! ' ' ! # & ! K # = X ; ( ) ) $ " C
( =
31
& ' ( $ ) *& * ) ' &' " ( ) # & # ' # ' ! # %*+ K # ( K ) ; ( =
( ) K* ;
43
* $ " ) ( @# , ! ! & & ! & ! - . ! C * X ) # > )
#
47
< =)# ? # /! ! & ) & " & "0 1 -)2 40 % !
! " ) " W = D ) W # )
# ) & U& C > ; # g $
59
* ! 5
U
( K( C"D\;K
63
]!
Monografie Komitetu Automatyki i Robotyki PAN
64
Nasze monografie
65
# ' ( #
statku morskiego Q 2 U2 - ( 4
VV
8
O
67
j < l XY2 - 0 8 % - # -8- Z[Y\
68
j < l ^^ 8 8 % . 88. Z[Z[
_[
]!
4 % ! ` - # . '
_Z
$
73
]!
& ( ( '(
74
? j # -( +
Z[Y\
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
* "D* ; K
Drodzy Czytelnicy, ' ) #)
naukowo-technicznego Pomiary # ) " $ & C )
)
M! = = # )&
) #
# ]
$ ) # M # )
Z =# M K* ) =
# =# M
[&
" #
) # ]
) #
) ) $ ) # &
"M ! # # =#
$ = )
) M
) & C ) ) #
Z$ ) [ M ! ) ! ' )
#
#
= ^ A ) ) #$ ZÂ ([ $ Z
DD[& C )
) A M
A
# # & ? $ !] ) #) ) !
= ^ # =
A OGUc OGOG& % ) ^ ^) M !
lB )M
#$ C E " + = D
] # ) M ]
# ] ! ) # & ? ) # ? A ( > J# A ) = K # #
# ) ) M K ? A = ] ) # =
numeru.
Redaktor naczelny kwartalnika Pomiary Automatyka Robotyka
& $& !& ; < '
3
4
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
Pomiary Automatyka Robotyka, ISSN 1427-9126, R. 23, Nr 1/2019, 5–10, DOI: 10.14313/PAR_231/5
. / $ % 046 # / % 8 9 ( 20 (
$ 3 ) # 8;9 8 % ;< =9 & > % ? & 8 & 4 # 4 @ A % & 8 B !& !=!E+ " <
The paper presents the implementation of the Fractional Order PID controller to control the temperature in the isothermal room applied in a pharmaceutical factory. The formula of the controller dedicated to transfer function model of the temperature is proposed, the stability analysis using the Matignon Theorem is also presented. Results of simulations show that the proposed controller is able to assure the better control quality than PID controller tuned with the use of auto-tuning function. 8 ` / $ % 046 & # % ( & F & $G8 H
1. Introduction The Fractional Order PID control (FOPID control) is the “classic� application of fractional order calculus in automation. Advantages of this idea are well known. A FOPID controller is able to assure better control quality and robustness than Integer Order (IO) PID due to use two additional parameters to its tune. FOPID controllers are considered by many Authors, for example: [1, 2, 8, 10, 13]. The recent state of the art is given in [8]. Tuning rules for FOPID controllers are given for example in [3, 7, 12, 14]. This paper is devoted to present the use of FOPID controller to temperature control in isothermal room, used for example in pharmaceutical factory. The job of the control system is to keep a constant, relatively high temperature inside the room. The main disturbation is variable outside temperature, varying dependently on season and weather. The paper is organized as follows: at the beginning some elementary ideas from fractional calculus are recalled. Next the considered control plant and its integer order transfer function model are given. Furthermore the FOPID controller is proposed and stability analysis for selected order of the controller is presented. Results are by simulations illustrated.
2. Preliminaries
,
Y2 (The elementary non integer order operator) The fractional order, integro-differential operator is defined as follows:
()
⎧ dÎą f t ⎪ ⎪ dt Îą ⎪ Îą a Dt f t = ⎨ f t ⎪ ⎪t ⎪ âˆŤ f Ď„ dĎ„ ⎊a
()
Îą >0
()
Îą =0
( )( )
Îą
(1)
Îą <0
where a and t denote time limits for operator calculation, a ∈ denotes the non-integer order of the operation. Next the Gamma Euler function is recalled (see for example [5]): ,
Z2 (The Gamma function)
()
∞
Γ x = âˆŤ t x −1e −tdt
(2)
0
Mittag-Leffler function is a non-integer order generalization of exponential function e lt and it plays crucial role in solution of FO state equation. The one parameter Mittag-Leffler function is defined as follows:
2.1. Elementary ideas A presentation of elementary ideas is started with a definition of a fractional order, integro-differential operator [2, 4, 5, 11]:
,
v2 (The one parameter Mittag-Leffler function)
. ' `
" # $ % & ' ( % . ( % )) !* )!*+ & % % *- ! )!*+
!
∞
k
( ) ∑ Γ kxÎą + 1 ( )
EÎą x =
(3)
k =0
The fractional-order, integro-differential operator is described by different definitions, given by GrĂźnwald and Letnikov (GL definition), Riemann and Liouville (RL definition) and Caputo (C definition). In the further consideration only C definition is applied:
5
The fractional order PID control of the forced air heating system ,
x2 (The Caputo definition of the FO operator)
C 0
()
DtÎą f t =
∞
1
âˆŤ Γ (N − Îą ) 0
Îą +1− N
dτ ,
(4)
⎛ ωh ⎞ âŽ&#x; âŽ? ωl âŽ
Ρ =⎜
where N – 1 < a < N denotes the fractional order of operation and (..) is the complete Gamma function expressed by (2). For the C definition the Laplace transform is defined as follows:
$ ⎥⎣ C0 Dtι
Z2v2 7 &+,
() ( ) Îą <0 f (t ) ⎤⎌ = s F (s ) − ∑ s n −1
Îą
Îą − k −1 0
()
Dtk f 0 ,
k =0
Îą > 0,
Consider the closed loop control system shown in the figure 1. It contains controller described by transfer function Gc(s) and plant described by the transfer function G(s). Denote the transfer function of the whole closed-loop control system by Gcl(s).
(5)
n −1 < Îą ≤ n ∈N
()
It equals to: Gcl s = The use of C definition (4) and (5) allows to define the fractional order transfer function, analogically to integer order case:
( ) = G (s )G (s ) . The FOPID conR ( s ) 1 + G ( s )G ( s )
Y s
c
c
troller is described by the following transfer function (see for example [3], page 278):
Gc(s) = kp + kis– + kds ,
bM s βM + ! + b1s β1 + b0 , aN s ιN + ! + a1s ι1 + a0
()
G1 s =
(10)
1 âˆ’Îą N
In (10) wl and wh describe the range of angular frequency, for which the approximation is expected to work properly. The steady state gain kf is set to assure the convergence the step response of the approximation to step response of the real plant in a steady state.
,
X2 (The Laplace transform for Caputo operator) $ ⎥⎣ C0 Dtι f t ⎤⎌ = s ι F s ,
Îą
⎛ ω ⎞N Îł =⎜ hâŽ&#x; âŽ? ωl âŽ
(N )
(Ď„ ) (t − Ď„ ) f
where:
(11)
(6) where kp, ki and kd are coefficients describing the proportional, integral and derivative actions respectively, a, b ∈ denote fractional orders of integral and derivative actions respectively.
where bM‌b1 and aN‌a1 are fractional orders, bM‌b0 and aN‌a0 are coefficients. In the particular case the commensurate order system can be considered. This denotes that: bm = mg, m = 1, ‌, M, aN = ng, n = 1, ‌, N, M ” N. Consequently the transfer function (6) takes the following form:
()
G1 s =
( ) (s )
bM s Îł
M
+ ! + b1s Îł + b0
Îł
N
+ ! + a1s Îł + a0
aN
( ). = D (s ) L sÎł
Fig. 1. Closed-loop control system Rys. 1. Układ zamknięty z ujemnym sprzęşeniem zwrotnym
(7)
Îł
Z2x2 % % The fundamental Bounded Input-Bounded Output (BIBO) stability condition for FO system described by the commesurate FO transfer function (7) is given by Matignon [6]. It is given also in [1 pages 21–22].
The transfer function (7) is quotient of quasi-polynomials L(s g) and D(s g). This makes easier the stability analysis. The transfer function (6) can be modeled at MATLAB platform only using its integer order, finite dimensional approximation. These approximations, both time-continuous and time-discrete are known. In this paper the Oustaloup Recursive Approximation (ORA) is applied. It is presented in the next subsection.
Theorem 1. (Matignon [6]) The fractional order system described by the transfer function (7) is stable if and only if:
2.2. The Oustaloup Recursive Approximation
( )
arg Îťn
The method proposed by [9] allows to approximate an elementary FO transfer function s a by the finite, integer-order transfer function, close to Pade approximation: N
sÎą ≈ k f âˆ? n =1
1+ 1+
s
Îźn s
=
νn
( ). (s )
ν n = ÎźnÎł , n = 1,! , N Îźn +1 = ν nΡ , n = 1,! , N − 1
O
M
I
A
R
Y
(12)
Z2X2 4 ' %
Ο1 = ωl Ρ
P
∀n = 1,! , N ,
2
It is important to note that for 0.0 < l < 1.0 the stability area covers also the positive part of complex plane and for a = 1.0 the stability area is the same, as for integer order system. The condition (12) can be easily examined using numerical methods. This will be shown in the “Simulations� section.
(8)
In (8) N is the order of the approximation, kf is the steady state gain, mn and nn are calculated as follows:
6
γπ
where ln is the n-th root of polynomial D(s g).
LORA s DORA
>
The control performance of the closed-loop control system can be tested using typical cost functions: settling time Ts and overshoot ymax expressed in % of a reference value. As the
(9)
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
integral cost function the Integral of Absolute Error (IAE) cost function will be used: IAE =
∞
âˆŤ e (t ) dt .
(13)
0
v2 4 ' v2Y2 4 % # The considered forced heating system is required to assure the constant, relatively high temperature in the room with cubature about 75 m3. The system contains only heater. The cooler is not necessary due to environmental conditions and required range of preset temperatures. The simplified scheme of the system is illustrated by the SCADA screen shown in the figure 2.
to avoid problems associated with the heterogeneous spatial temperature distribution in the room. Trends of temperature and control caused by the change of the preset temperature from 28 °C to 37 °C are shown in the figure 3. These trends were measured during normal work of the heating system, controlled by classic PID controller. This “passive� method was only possible to apply in the considered case. The sample time during experiments was equal: h = 30 s. The outside temperature is the biggest disturbance for the system. However, in case of such a short period of time it can be assumed as constant.
v2Z2
4 % %
% '
in the room The second order IO transfer function model of the temperature was obtained using the data shown in the figure 3 and tfest MATLAB function. It is as follows:
()
G s =
0.0002493s + 2.117e −07 s + 0.002379s + 5.053e −07 2
(14)
The fitting the response of model to response of the real plant is shown in the figure 4.
Fig. 2. The forced air heating system Rys. 2. System ogrzewania powietrza
The cold air from the room is sucked via return duct to the heater and next, after being heated by the heater 2E1 and filtered by the HEPA filter, is supplied to the room via supply duct. The filter pollution is monitored by the pressure switch 11PS2. The electric heater 2E1 is protected from damage by monitoring air flow with the use of another pressure switch. The control signal is the power of the heater 2E1 expressed in the percent of its maximal range. The rotational speed of the blower 2M1 can be also applied as a control signal, but in this case, it is maximal and constant. The controlled temperature is measured inside the return duct by the sensor Pt 1000 (element 6TE1 in the scheme). This location of the sensor allows
Fig. 4. The comparison the model (14) to real data Rys. 4. PorĂłwnanie modelu (14) z danymi rzeczywistymi
x2 The general form of the transfer function of the plant with respect to (14) is expressed as follows:
()
G s =
b1s + b0 . s 2 + a1s + a0
(15)
If we recall the transfer function of the controller (11), then the open-loop transfer function Go(s) of the control system is as follows:
()
Go s =
where:
Fig. 3. The temperature measured by sensor 6TE1 and the power of heater 2E1 Rys. 3. Temperatura zmierzona przez czujnik 6TE1 i moc podgrzewacza 2E1
( ) =ds M (s ) L s
5
ι + β +1
+ d 4s ι +1 + d 3s ι + β + d 2s ι + d1s + d 0 , (16) s ι + 2 + a1s ι +1 + a0s ι
⎧d 5 ⎪ ⎪d 4 ⎪ ⎪d ⎪ 3 ⎨ ⎪d 2 ⎪ ⎪d1 ⎪ ⎪⎊d 0
= b1kd = b1k p = b0kd = bok p
(17)
= b1ki = boki
7
The fractional order PID control of the forced air heating system The closed-loop transfer function of the control system is following:
()
Gcl s =
()
( ), L (s ) + M (s ) D (s ) L s
=
L s
that the values of coefficients kp, ki, kd for the FOPID are much bigger, than suitable values for classic PID. This is caused by the fact that the steady-state gain of the ORA approximation applied to describe FO elements in the FOPID is very small. This causes that the work of the FO parts must be amplified by bigger values of coefficients. The SIMULINK model contains also the control system with typical PID controller, tuned using the available auto-tune function. It is employed to compare the proposed results to typical, known solution. The work of the system during change of preset temperature from 27 ºC to 37 ºC was tested. The level of the control signal in the both control systems was limited from 0% to 100%. Parameters of the both controllers are collected in the table 2. Results are illustrated by the figure 6, the performance indices: settling time, overshoot and IAE (13) are given in the table 3. The roots of characteristic quasi-polynomial for FOPID described in the table 2 are shown in figure 7.
(18)
where L(s) and M(s) are given by (16), D(s) is the characteristic polynomial of the closed-loop control system:
(19) Next assume that orders of the both fractional order parts are equal: a = b = 0.75. It allows to substitute: p = s 0.25 and simplifies the further analysis. Additionally results of simulations show that this value of a and b assures very good control performance. Consequently the characteristic quasi-polynomial D(p) = D(s0.25) takes the following form:
Table 2. Parameters of controllers Tabela 2. Parametry regulatorów
(20) where d5 – d0 are given by (17). Consequently the stability condition (12) turns to the following form:
( )
arg λn
π
>
∀n = 1,! ,11,
8
Controller
Parameters
PID
kp = 26.9920, ki = 0.0772, kd = 289.6436, N = 0.9142
FOPID
kp = 500.0, ki = 0.1, kd = 200.0, a = 0.75, b = 0.75
(21)
where ln is the n-th root of the quasi polynomial (20). The roots of quasi polynomial (20) can be calculated numerically for given controller parameters kp, ki and kd. This will be shown in the next section.
X2 Simulations were executed with the use of transfer function model of plant (14) and controller (11). The fractional order parts of the controller were modeled using ORA approximation (8). The parameters of ORA approximation applied in experiments are given in the table 1. The SIMULINK model Table 1. The ORA approximation parameters Tabela 1. Parametry aproksymacji ORA
N
wl
wh
8
10–4
104
Fig. 6. PID vs FOPID controller step responses Rys. 6. Porównanie odpowiedzi skokowych układu z regulatorem PID i regulatorem FOPID
Fig. 5. The SIMULINK model applied to tests Rys. 5. Modele simulinkowe zastosowane do testów Fig. 7. The roots location for quasi polynomial (20) and parameters of FOPID from table 2 in the complex plane Rys. 7. Lokalizacja na płaszczyźnie zespolonej pierwiastków quasi wielomianu (20) dla układu z regulatorem FOPID o parametrach zamieszczonych w tabeli 2
applied to tests is shown in the figure 5. The FOPID was tuned experimentally to obtain the shortest possible settling time with the smallest overshoot. In the table 1 it can be noted
8
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
Table 3. Performance indices for the both control systems Tabela 3. Wskaźniki jakości regulacji dla obu systemów sterowania Controller
Settling time [s]
Overshoot [% of
IAE (13)
preset value]
PID
734.82
0.0
2.9304e+03
FOPID
540.37
3.05
2.8922e+03
V2 The main conclusion from the paper is that the proposed simple FOPID controller is able to assure the better control performance than typical PID controller tuned with the use of the auto-tune function. The FOPID controller is easy to tune using simple transfer function model of the plant and simulations. The further investigations will cover a synthesis of robust FOPID controller for the considered plant. Another interesting problem is to give general stability conditions for the considered FOPID control system. The practical verification at the real control system is recently impossible due to running production in the factory. Experiments will be run as soon as possible, probably during longer service break.
! % 1. Caponetto R., Dongola G., Fortuna L., Petras I., Fractional order systems: Modeling and Control Applications. [in:] Chua L.O., editor, World Scientific Series on Nonlinear Science, 1–178. University of California, Berkeley, 2010. 2. Das S., Functional Fractional Calculus for System Identification and Control. Springer, Berlin 2010. 3. Faieghi M.R., Nemati A., On fractional-order PID design. [in:] Michalowski T., editor, Applications of MATLAB in Science and Engineering, 273–292. InTech, Rijeka Croatia, 2011. 4. Kaczorek T., Singular fractional linear systems and electrical circuits. “International Journal of Applied Mathematics and Computer Science�, Vol. 21, No. 2, 2011, 379–384. 5. Kaczorek T., Rogowski K., Fractional Linear Systems and Electrical Circuits. Bialystok University of Technology, Bialystok, 2014. 6. Matignon D., Stability results for fractional differential equations with applications to control processing. [in:] IMACS-SMC Proceedings, Lille, France, July 1996, 963– 968.
7. Merrikh-Bayat F., Mirebrahimi N., Khalili M.R., Discrete-time fractional-order PID controller: Definition, tuning, digital realization and some applications. “International Journal of Control, Automation, and Systemsâ€?, Vol. 13, No. 1, 2015, 81–90. 8. Edet E., Katebi R., On Fractional-Order PID Controllers, Fractional order systems: Modeling and Control Applications. [in:] Preprints of the 3rd IFAC Conference on Advances in Proportional-Integral-Derivative Control, Ghent, Belgium, May 9–11, 2018, 739–744. 9. Oustaloup A., Levron F., Mathieu B., Nanot F., Frequency-band complex noninteger differentiator: characterization and synthesis. IEEE Transactions on Circuits and Systems I: Fundamental Theory Applications, Vol. 47, No. 1, 2000, 25–39, DOI: 10.1109/81.817385. 10. Petras I., Fractional – Order feedback control of a DC motor. “Journal of Electrical Engineeringâ€?, Vol. 60, No. 3, 2009, 117–128. 11. Podlubny I., Fractional Differential Equations. Academic Press, San Diego 1999. 12. Ranganayakulu R., Babu G.U.B., Rao A.S., Patle D.S., A comparative study of fractional order PI /PI DÎź tuning rules for stable first order plus time delay processes. “Resource-Efficient Technologiesâ€?, Vol. 2, No. 1, 2016, 136–152. 13. Teplakov A., Fractional-order Modeling and Control of Dynamic Systems. PhD thesis, Tallin Unviersity of Technology, Estonia, Tallin 2015. 14. Valerio D., Costa J., Tuning rules for fractional PID controllers. [in:] Proceedings of the 2nd IFAC Workshop on Fractional Differentiation and its Applications Porto, Portugal, July 19–21, 2006, 1–6.
046 ( % % ( 4 Artykuł prezentuje implementację regulator PID ułamkowego rzędu (FOPID) do sterowania temperaturą w pomieszczeniu izotermicznym stosowanym w fabryce farmaceutycznej. Zaproponowano formułę regulatora dla modelu obiektu opisanego transmitancją oraz analizę stabilności z wykorzystaniem Tw. Matignona. Wyniki badań symulacyjnych wskazują, şe proponowany regulator zapewnia lepszą jakość regulacji, niş typowy regulator PID dostrojony zuşyciem autotuningu. ' ` ( 046& ( & F I& $G8
9
The fractional order PID control of the forced air heating system
20 (
$ 3 !5 54
) # 3 4 6/
' ( %
% % E'(
9 F " *+J- 9 F % B *+KK& 0 6 % 6 8 % G F *++E % )!!+ 8;9 # % . = ( L" & 0 %M 9 F ( 8;9 6 # 8 = *+KK& # 4 )!*)N )!*J % % # # # ? = & 8 & 4 # % A = % ? ( ( 8;9 9 O % & # % % ( % = & & % & 0P % 868 = & F F
9 F *Q !J *++- . B # % ( A # ? ( ( % ( % % R )!*K / # ? ? ( (& 8 = & % A % ? ( = ( 8;9 # % . ( 9 0 6 % 6 # 8 % G F = 9 ( # % = 9 O % ( # % ( # % % # %
10
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
Pomiary Automatyka Robotyka, ISSN 1427-9126, R. 23, Nr 1/2019, 11–16, DOI: 10.14313/PAR_231/11
> O S < % Jacek Michalski 1, Piotr Kozierski 2, 1 1
0 0 T & > % ? & 4 8 & G F 4 @ 4 # & % 8 G F & 0 & J!=+JE 0 T
2
0 0 T & > % 4 # & 4 8 G F & % %< ? 0 ( < & 0 & J!=+JE 0 T
Streszczenie: W artykule poruszono problem estymacji stanu systemów dynamicznych oraz zaproponowano nową metodę jego rozwiązania – wielokrotny filtr cząsteczkowy. Jest to odmiana filtru cząsteczkowego pozwalająca na zrównoleglenie jego pracy przez podział na niezaleşne filtry tak, by umoşliwić implementację algorytmu, takşe na urządzeniach o nieduşej mocy obliczeniowej. Algorytm został zaimplementowany dla obiektu jedno- oraz wielowymiarowego, a jakość estymacji porównano dla róşnej liczby cząsteczek. Do oceny działania algorytmu wykorzystano wskaźnik jakości aRMSE. Na podstawie badań stwierdzono, iş zrównoleglenie pracy filtru cząsteczkowego moşe poprawić działanie algorytmu. ( ` & % & O A & O S & ( A & O S
1. Wprowadzenie Estymacja stanu jest działem nauki bardzo często wykorzystywanym w praktyce, szczególnie w środowiskach pomiarowych, w których oddziałują czynniki stochastyczne. Umoşliwia ona zmniejszenie negatywnego wpływu oddziaływania szumów pomiarowych, szczególnie gdy liczba sygnałów pomiarowych przekracza liczbę zmiennych stanu. Estymacja jest szeroko stosowana w wielu dziedzinach techniki, na przykład w analizie ruchu robotów latających UAV (ang. Unnamed Aerial Vehicle) [1, 2], robotów mobilnych [3], czy teş na większą skalę, w analizie działania sieci elektroenergetycznej [4]. W ostatnich latach bardzo waşną rolę w estymacji stanu odgrywa filtr cząsteczkowy PF (ang. Particle Filter), którego pierwszy algorytm BF (ang. Bootstrap Filter) zaproponowali w 1993 r. Gordon, Salmond i Smith [5]. Filtr cząsteczkowy coraz częściej wypiera algorytmy Kalmana KF (ang. Kalman Filter) [6] ze względu na lepszą jakość estymacji dla systemów nieliniowych [7]. Algorytmy PF, jako jedne z niewielu metod, radzą sobie z estymacją szumów o dowolnej funkcji gęstości prawdopodobieństwa PDF (ang. Probability Density Function),
podczas gdy algorytmy KF przeznaczone są do estymacji jedynie systemów z zakłóceniami o rozkładzie Gaussa. Przez lata powstała duşa liczba modyfikacji filtrów cząsteczkowych, np. inteligentny filtr cząsteczkowy IPF (ang. Intelligent Particle Filter) [8], czy iteracyjny rozszerzony cząsteczkowy filtr Kalmana IEKPF (ang. Iterated Extended Kalman Particle Filter) [9]. Cały czas tworzone są nowe algorytmy estymacji stanu, jak na przykład metoda oparta na sztucznych sieciach neuronowych [10] lub systemach rozmytych [11]. W niniejszej pracy autorzy skupili się na filtrze cząsteczkowym ze względu na jego bardzo dobre właściwości w odniesieniu do estymacji stanu obiektów nieliniowych. Zaproponowana została nowa metoda – wielokrotny filtr cząsteczkowy MultiPDF PF (ang. Multi Probability Density Function Particle Filter) oraz przedstawione zostało jego działanie dla róşnego typu obiektów nieliniowych. W drugim rozdziale sformułowano rozwiązywany problem, w trzecim rozdziale przedstawiono filtr cząsteczkowy oraz zaproponowano modyfikacje. W czwartym rozdziale pokazano modele badanych obiektów. Rozdział piaty prezentuje wykresy obrazujące wyniki symulacji. W rozdziale szóstym zawarto wnioski końcowe.
Z2 % ( ' . ' `
R B & +E' . ( % )E !) )!*+ & % % )+ ! )!*+
!
Rozwaşany jest układ dynamiczny, nieliniowy, dyskretny, dany równaniami stanu
(
)
⎧x (k +1) = f x (k ) , u (k ) ; k + v (k ) ⎪⎪ ⎨ (k ) (k ) (k ) ⎪ y =h x +n ⎪⎊
( )
(1)
11
? E ] ) # )M ) gdzie x(k) jest wektorem stanu w k-tym kroku czasowym, u(k) – wektorem wejść, y(k) – wektorem pomiarów (wyjść), v(k) to wektor szumów wewnętrznych, a n(k) to wektor szumów pomiarowych; f jest wektorem funkcji przejść, a h – wektorem funkcji pomiarowych.
rytmu SIR staje się równieş bardziej skomplikowane. Algorytm Bootstrap przedstawia się następująco:
Wielokrotny filtr cząsteczkowy ma za zadanie w kaşdym kroku estymować wartości wektora stanu x(k), które nie muszą być mierzalne, mając dane sygnały wejściowe i pomiarowe. Na podstawie estymowanego wektora stanu oblicza się estymowane wartości pomiarów, których wartości powinny być przefiltrowane (ze zmniejszonym wpływem szumów pomiarowych).
2. Prognoza. Wylosowanie czÄ…steczek z modelu przejĹ›cia xi,(k) ~ p(x(k)|xi,(k−1)).
1. Inicjalizacja. Wylosowanie Np wartości początkowych xi,(0) z początkowego PDF p(x(0)), ustawienie numeru iteracji k = 1.
3. Aktualizacja. Obliczenie wagi kaĹźdej z nowych czÄ…steczek na podstawie modelu pomiarowego
(
i ,(k ) (k ) i ,(k ) q � p y |x
v2 Filtry cząsteczkowe są przeznaczone do estymacji stanu obiektów nieliniowych o szumach z dowolnym rozkładem gęstości prawdopodobieństwa, co jest ich zdecydowaną przewagą nad większością metod estymacji (jak np. filtry Kalmana) [7]. Podstawą działania kaşdego z algorytmów PF jest rekursywny filtr Bayesa dany wzorem (2)
(
p x
|Y
(k )
)= (
p y
(k )
) â‹… p (x ( ) | Y ( ) ) , () ( ) p (y | Y )
|x
(k )
k −1
k
k
k
k −1
k
(
p x
(k )
|Y
(k )
)
)
(
(k ) (k ) = pˆ x | Y =
Np
∑q
()
i, k
(3)
)
(
I
A
R
Y
•
A
U
T
(6)
j ,(k ) q
()
2. Obliczanie wartoĹ›ci estymat w kaĹźdym b-tym (b = 1, ..., Nf) z rĂłwnolegle dziaĹ‚ajÄ…cych filtrĂłw. − Prognoza. Wylosowanie czÄ…steczek z modelu przejĹ›cia
(4)
()
O
M
(
p x
(k )
(
)
i , k −1
| xb
).
− wagi kaĹźdej z nowych czÄ…steczek − Aktualizacja. Obliczenie na podstawie modelu pomiarowego
PoniĹźej przedstawiono algorytm Bootstrap filtru czÄ…steczkowego [5], ktĂłry stanowiĹ‚ podstawÄ™ dla zaproponowanej w pracy metody wielokrotnego filtru czÄ…steczkowego. Jest on pewnym uproszczeniem w stosunku do podstawowego algorytmu filtru czÄ…steczkowego – sekwencyjnej funkcji waĹźnoĹ›ci z powtĂłrnym prĂłbkowaniem SIR (ang. Sequential Importance Resampling) [13], ktĂłry zakĹ‚ada losowanie czÄ…steczek w etapie prognozy nie z modelu przejĹ›cia, a z funkcji waĹźnoĹ›ci g(x(k)|xi,(k−1), y(k)). W zwiÄ…zku z tym obliczanie wag czÄ…steczek za pomocÄ… algo-
M
j =1
i, 0
gdzie d(¡) oznacza deltę Diraca.
O
Np
1. Inicjalizacja. Wylosowanie Np/Nf wartości początkowych xb z początkowego PDF p(x(0)), ustawienie numeru iteracji k = 1.
i, k
P
∑
Metoda wielokrotnego filtru cząsteczkowego jest to zaproponowana modyfikacja podstawowego algorytmu. W myśl MultiPDF PF w kaşdym kroku symulacji działa równolegle i niezaleşnie od siebie Nf filtrów cząsteczkowych obliczających estymaty dla tego samego obiektu, z których kaşdy ma odpowiednio mniejszą liczbę cząsteczek – Np/Nf, gdzie Np jest liczbą cząsteczek całego algorytmu. Na końcu kaşdego kroku wyniki z poszczególnych Nf filtrów są uśredniane i w ten sposób obliczana jest estymata w k-tym kroku próbkowania. Dzięki podziałowi algorytmu na filtry o odpowiednio mniejszej liczbie cząsteczek, metodę MultiPDF PF moşna implementować na urządzeniach niedysponujących odpowiednio duşą mocą obliczeniową do zaimplementowania tradycyjnych metod filtru cząsteczkowego. Dla parametru Nf = 1 algorytm MultiPDF Particle Filter staje się tradycyjnym filtrem cząsteczkowym Bootstrap. Ponişej przedstawiono algorytm wielokrotnego filtru cząsteczkowego. Algorytm ten naleşy powtórzyć dla kaşdego z Nf równolegle działających filtrów cząsteczkowych.
xb
12
i ,(k ) q
v2Z2 %
k
k i, k δ x( ) − x ( ) ,
i =1
=
Moşna stworzyć dowolną liczbę filtrów cząsteczkowych przez odpowiedni dobór funkcji waşności. Znane są róşne rodzaje filtrów cząsteczkowych [13–16].
gdzie Y(k) jest zbiorem pomiarĂłw wyjĹ›cia z k poczÄ…tkowych krokĂłw, p(x(k)|Y(k)) jest to funkcja gÄ™stoĹ›ci prawdopodobieĹ„stwa a posteriori, p(y(k)|x(k)) oznacza wiarygodność, czyli prawdopodobieĹ„stwo warunkowe wystÄ…pienia pomiaru y(k) pod warunkiem, Ĺźe wektor stanu przyjmuje wartość x(k), p(x(k)|Y(k-1)) to rozkĹ‚ad prawdopodobieĹ„stwa a priori, a p(y(k)|Y(k-1)) stanowi współczynnik normujÄ…cy [12]. Autorzy do obliczania wag zastosowali wzĂłr bez współczynnika normujÄ…cego, poniewaĹź w odpowiednim, późniejszym kroku algorytmu wagi zostajÄ… znormalizowane tak, by ich suma byĹ‚a rĂłwna jednoĹ›ci. W kaĹźdym algorytmie filtru czÄ…steczkowego PDF jest reprezentowana przez zestaw czÄ…steczek, z ktĂłrych kaĹźda i-ta czÄ…steczka jest zĹ‚oĹźona z wektora stanu xi oraz wagi qi. Waga czÄ…steczki stanowi prawdopodobieĹ„stwo przyjmowania przez niÄ… wartoĹ›ci prawdziwej wektora stanu, ktĂłra to wartość bardzo czÄ™sto w praktyce jest niemierzalna. KaĹźda i-ta czÄ…steczka moĹźe zostać zapisana jako para {xi, qi}. Gdy liczba czÄ…steczek Np jest wystarczajÄ…co duĹźa, PDF a posteriori moĹźna zapisać N p →∞
()
i, k
5. Resampling systematyczny (opisany na końcu rozdziału). 6. Koniec iteracji. Obliczenie estymaty x(k), aktualizacja numeru iteracji k := k+1, przejście do kroku 2.
(2)
k −1
k
q
k −1
k
( ( ) | Y ( ) ) = âˆŤ p ( y ( ) | x ( ) ) â‹… p ( x ( ) | Y ( ) ) dx ( ) ,
p y
(5)
4. Normalizacja. Skalowanie wag w taki sposób, by ich suma była równa 1.
v2Y2 . # ~ '
(k )
)
(
i ,(k ) i ,(k ) (k ) q b � p y | xb
)
(7)
− Normalizacja. Skalowanie wag w taki sposĂłb, by ich suma byĹ‚a rĂłwna 1. ()
i, k
qb
A
T
Y
K
A
•
R
O
B
()
i, k
=
O
q b
∑
Np
(8)
() q
j, k j =1 b
T
Y
K
A
NR 1 / 20 1 9
− Resampling systematyczny (opisany na koĹ„cu rozdziaĹ‚u). (k ) − Koniec iteracji. Obliczenie estymaty xb . 3. Obliczenie estymaty x(k) jako Ĺ›redniej arytmetycznej Nf wynikĂłw uzyskanych w poprzednim kroku, aktualizacja numeru iteracji k := k + 1, przejĹ›cie do kroku 2.
noległe części. Modele badanych obiektów zostały przedstawione ponişej. • Ob1
v2v2 ! ' # Resampling polega na ponownym wylosowaniu wartoĹ›ci czÄ…steczek, ale tylko spoĹ›rĂłd tych juĹź istniejÄ…cych. PrawdopodobieĹ„stwo wylosowania danej czÄ…steczki jest rĂłwne wartoĹ›ci jej znormalizowanej wagi. Po resamplingu wagi wszystkich czÄ…steczek sÄ… sobie rĂłwne i wynoszÄ… q i = N p−1 (a w przypadku filtru wielokrotnego qi = Nf/Np). PoniĹźej przedstawiono algorytm resamplingu systematycznego [17] wykorzystany w aktualnych badaniach. Dany algorytm powtarza siÄ™ w kaĹźdym kroku symulacji. Dla wielokrotnego PF naleĹźy w kaĹźdym miejscu wartość Np zamienić na Np/Nf. WiÄ™cej o resamplingu i róşnych jego rodzajach moĹźna znaleźć w [17].
(k ) ⎧ (k +1) 25x (k ) (k ) ⎪x = 0,5x + + 8 cos 1,2k + v 2 (k ) ⎪⎪ 1+x ⎨ 2 ⎪ (k ) (k ) x (k ) ⎪ +n y = ⎪⎊ 20
(
v
(k )
(
)
& 0; 10 , x
() 0
n
(k )
)
( )
& 0; 1
= 0,1
• Ob2
1. Ustawienie j = 1, i = 1, Q = qi,(k). Wylosowanie u ~ -[0, 1]/Np. 2. Dopóki Q < u następuje inkrementacja j = j + 1, Q = Q + qj,(k). () 3. Zapamiętanie cząsteczki xj,(k) jako wylosowanej x . Ustai,(k) wienie wagi otrzymanej cząsteczki jako q = 1/Np. i, k
4. Inkrementacja wartości i = i + 1 oraz u = u + 1/Np. Przejście do kroku 2.
x2 ~
(k )
Dla uproszczenia zapisu załoşono, şe x i Do badań wykorzystano dwa obiekty, z czego oba są autonomiczne (nie mają sygnałów wejściowych). Pierwszy z nich, Ob1, jest bardzo często wykorzystywany do badania filtrów cząsteczkowych [13, 18]. Ze względu na swoją silną nieliniowość bardzo dobrze pokazuje przewagę filtrów cząsteczkowych nad algorytmami wykorzystującymi linearyzację punktową [7]. Drugi system Ob2 jest obiektem wielowymiarowym, zaproponowanym przez autorów na potrzeby badań nad MultiPDF PF, poniewaş pokazuje on korzyści płynące z zastosowania podziału algorytmu na pracujące niezaleşnie od siebie, rów-
(k )
( ) & ( 0; 4 ) ,
v1:5 & 0; 1 , (k )
n6:10
x
(0)
(k )
(k )
≥ xi + 5 .
( ) & ( 0; 1) ,
n1:6 & 0; 6 , (k )
n11:15
T
= ⎥⎣0 0 0 0 0 ⎤⎌ .
Zastosowane oznaczenia: (m, s2) – rozkład normalny o wartości oczekiwanej m i wariancji s2.
Rys. 1. Wartości aRMSE dla Ob1 Fig. 1. Values of aRMSE for Ob1
13
? E ] ) # )M )
aRMSE dla Ob2 Fig. 2. Values of aRMSE for Ob2
X2 <
zaproponowany algorytm polepsza jakość estymacji działając z nieliniowymi obiektami wielowymiarowymi. Warto zaznaczyć, şe wąskie gardło w algorytmach PF stanowi obliczanie wag oraz resampling, są tam bowiem potrzebne informacje o wartościach wszystkich cząsteczek. Dzięki niezaleşnemu działaniu kaşdego z filtrów równoległych, kaşdy z nich ma niezbędne informacje, bez konieczności wzajemnego komunikowania się. Dzięki temu zaproponowany algorytm ma jednocześnie potencjał do przyspieszenia działania w porównaniu do algorytmu standardowego. Jak wynika z przeprowadzonych badań, oprócz moşliwości implementacji na słabszych obliczeniowo urządzeniach, zaproponowany algorytm moşe takşe w odpowiednich przypadkach poprawić jakość estymacji. Autorzy w przyszłości planują przetestować działanie wielokrotnego filtru cząsteczkowego dla róşnego typu obiektów (równieş posiadających sygnały wejściowe – sterujące) aby ocenić, w jakich przypadkach algorytm poprawia jakość estymacji. Planuje się takşe przeprowadzić badania nad modyfikacją wewnętrznej struktury algorytmu MultiPDF, na przykład zmieniając funkcję waşności, czy teş sposób obliczania estymaty na podstawie wyników z poszczególnych, równoległych filtrów.
Symulacje przeprowadzono kaşdorazowo dla M = 1000 kroków czasowych. Kaşda symulacja Ob1 została powtórzona co najmniej 1000 razy, dla Ob2 z kolei co najmniej 200 razy. Za kaşdym razem przebiegi wzorcowe generowane były na nowo, w sposób losowy. Do oceny jakości estymacji został wykorzystany wskaźnik jakości aRMSE dany ponişszymi wzorami.
(9)
RMSE i = MSE i
aRMSE =
1 Nx
(10)
Nx
∑ RMSE
(11)
i
i =1
gdzie M to długość (ksymulacji (liczba kroków), Nx to liczba ) to wartości i-tej zmiennej stanu zmiennych stanu, xˆi i w k-tym kroku próbkowania, odpowiednio estymowanej oraz prawdziwej. Ponişsze wykresy przedstawiają wartości wskaźników jakości z 95% przedziałem ufności (naniesione zostały odchylenia standardowe 2s w kaşdą stronę) [19].
~ # % 1. Hajiyev C., Soken H.E., Robust adaptive Kalman filter for estimation of UAV dynamics in the presence of sensor/ actuator faults, “Aerospace Science and Technology�, Vol. 28, No. 1, 2013, 376–383, DOI: 10.1016/j.ast.2012.12.003. 2. Marantos P., Koveos Y., Kyriakopoulos K.J., UAV State Estimation using Adaptive Complementary Filters, “IEEE Transactions on Control Systems Technology�, Vol. 24, No. 4, 2016, 1214–1226, DOI: 10.1109/TCST.2015.2480012. 3. Chang C., Ansari R., Khokhar A., Multiple object tracking with kernel particle filter, IEEE Computer Society Conference on Computer Vision and Pattern Recognition, June 2005, Vol. 1, 566–573, DOI: 10.1109/CVPR.2005.243. 4. Korres G.N., A Distributed Multiarea State Estimation, “IEEE Transactions on Power Systems�, Vol. 26, No. 1, 2011, 73–84, DOI: 10.1109/TPWRS.2010.2047030. 5. Gordon N.J., Salmond D.J., Smith A.F.M., Novel Approach to Nonlinear/non-Gaussian Bayesian State Estimation,
6. Wnioski Zbadano działanie zaproponowanego nowego algorytmu dla róşnych obiektów, przy podziale na róşną liczbę równoległych filtrów. Bazując na wynikach symulacji moşna stwierdzić, iş dla niektórych typów obiektów podział na odpowiednią, niewielką liczbę równoległych filtrów poprawia jakość działania algorytmu. W przypadku Ob1 jakość estymacji pogarszała się wraz ze wzrostem parametru Nf. Natomiast dla Ob2 przy podziale na 2 i 3 równoległe filtry została zaobserwowana poprawa wartości wskaźnika aRMSE. Przy dalszych podziałach jakość estymacji ulegała pogorszeniu. Moşna na tej podstawie stwierdzić, şe
14
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
“IEE Proceedings F – Radar and Signal Processingâ€?, Vol. 140, No. 2, 1993, 107–113. DOI: 10.1049/ip-f-2.1993.0015. 6. Kalman R.E., A New Approach to Linear Filtering and Prediction Problems, “Journal of Basic Engineeringâ€?, Vol. 82, No. 35, 1960, 35–45. DOI: 10.1155/1.3662552. 7. Michalski J., Kozierski P., ZiÄ™tkiewicz J., PorĂłwnanie metod estymacji stanu systemĂłw dynamicznych, „Pomiary Automatyka Robotykaâ€?, Vol. 21, Nr 4, 2017, 41–47, DOI: 10.14313/PAR_226/41. 8. Yin S., Zhu X., Intelligent particle filter and its application to fault detection of nonlinear system, “IEEE Transactions on Industrial Electronicsâ€?, Vol. 62, No. 6, 2015, 3852–3861, DOI: 10.1109/TIE.2015.2399396. 9. Liang-Qun L., Hong-Bing J., Jun-Hu L., The iterated extended Kalman particle filter, [in:] IEEE International Symposium on Communications and Information Technology, ISCIT 2005, Vol. 2, 1213–1216, DOI: 10.1109/ISCIT.2005.1567087. 10. Yadaiah N., Sowmoya G., Neural Network Based State Estimation of Dynamical Systems, [in:] International Joint Conference on Neural Networks, IJCNN’06, 1042–1049, DOI: 10.1109/IJCNN.2006.246803. 11. Miranda V., Pereira J., Saraiva J.T., Load allocation in DMS with a fuzzy state estimator, “IEEE Transactions on Power Systemsâ€?, Vol. 15, No. 2, 2000, 529–534, DOI: 10.1109/59.867136.
12. Straka O., Simandl M., Particle Filter with Adaptive Sample Size, “Kybernetikaâ€?, Vol. 47, No. 3, 2011, 385–400. 13. Michalski J., Kozierski P., ZiÄ™tkiewicz J., Comparison of auxiliary and Likelihood Particle Filters for state estimation of dynamical systems, “PrzeglÄ…d Elektrotechnicznyâ€?, Vol. 94, 2018, 86–90, DOI: 10.15199/48.2018.12.19. 14. Arulampalam S., Maskell S., Gordon N., Clapp T., A Tutorial on Particle Filters for On-line Nonlinear/ Non-Gaussian Bayesian Tracking, “IEEE Transactions on Signal Processingâ€?, Vol. 50, No. 2, 2002, 174–188, DOI: 10.1109/78.978374. 15. Michalski J., Kozierski P., ZiÄ™tkiewicz J., Comparison of Particle Filter and Extended Kalman Particle Filter, “Studia z Automatyki i Informatykiâ€?, Vol. 42, 2017, 43–51. 16. Michalski J., Kozierski P., ZiÄ™tkiewicz J., Giernacki W., Likelihood Particle Filter and Its Proposed Modifications, “Studia z Automatyki i Informatykiâ€?, Vol. 43, 2018, 81–93. 17. Kozierski P., Lis M., ZiÄ™tkiewicz J., Resampling in Particle Filtering – Comparison, “Studia z Automatyki i Informatykiâ€?, Vol. 38, 2013, 35–64. 18. Kitagawa G., Monte Carlo Filter and Smoother for Non-Gaussian Nonlinear State Space Models, “Journal of Computational and Graphical Statisticsâ€?, Vol. 5, No. 1, 1996, 1–25, DOI: 10.2307/1390750. 19. Florek A., Mazurkiewicz P., SygnaĹ‚y i systemy dynamiczne, Wyd. 2, PoznaĹ„ 2015, ISBN: 978-83-7775360-6.
B 06/ 0 / # ? # 6 Abstract: In this paper the problem of state estimation of dynamical systems has been discussed and the new solution, named MultiPDF Particle Filter has been proposed. It is a modification of Particle Filter that allows to parallelize its work by dividing into independent filters in a way to enable the implementation of the algorithm also on devices with low computing power. The algorithm has been implemented for a one- and multi-dimensional object, and the quality of the estimation has been compared for a different number of particles. The quality index aRMSE has been used to evaluate the algorithm’s performance. Based on the simulation results it was found that the work parallelization of a Particle Filter can improve estimation quality of the algorithm. Keywords` & % & A / & 0 / & A / & B 06/ 0 /
" ; !
"
+E'
'(
% % ( ( 44 0 = 0 T F = S F < L M& % = % O " O S = ( & < @ S (
6 0 = 0 T ; < F = T % F < L % # S ( F I ( M = S( L & M
15
NR 3/2015
16
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
Pomiary Automatyka Robotyka, ISSN 1427-9126, R. 23, Nr 1/2019, 17–20, DOI: 10.14313/PAR_231/17
F U ? # / G V % ? ( F X F # G # 0 2 $ <= $ > # . ( & B / & 8 % A F K& !)=E)E > & 0 %&
# # Y "8 4?>4 G X N 4 % G 4 # 8 % B 0480& 8 R )!)& !)=-KJ > & 0 %
Abstract: The paper presents a new concept of video presentation for subjective quality evaluation methods with a reference given in parallel. The idea is to split a full resolution video into an n-picture matrix and encode each cell differently to compare aspects of processing such as encoding parameters or lossy compression algorithms. Conducted experiments show that it is possible to get more accurate results, shortening the time of evaluation at least by half with a less complicated and cheaper experimental station. 8 ` % & F Z
1. Introduction Subjective user-perceived video quality ratings provided by a human audience and expressed as a mean opinion score (MOS) are still considered to be the most reliable measure of compressed video quality [1, 2]. Aiming at unification, in its recommendations the International Telecommunication Union (ITU) describes the methods of conducting subjective experiments. There are methods that use single stimuli and others which use more than one picture to score quality. Additional pictures are used as a reference, which enables comparison and helps to formulate assessment on a given scale. When stimuli are displayed sequentially, such as in the case of pair comparison (PC) [3], degradation category rating (DCR) [3], or subjective assessment methodology for video quality (SAMVIQ) [4], the experiment takes more time and the measurement accuracy is decreased by memory effects [7]. This problem does not occur in the case of methods that use stimuli given in parallel, such as in the case of the simultaneous double stimulus for continuous evaluation (SDSCE) [5] or triple stimulus continuous evaluation scale (TSCES) [6], but in turn they require duplication of video encoders and monitors, which multiplies the price of the experimental station. This is a big drawback, especially in the case of new technologies such as 8K ultra-high-definition (UHD) television, for example. The other problem is synchronizing video signals and the unnatural viewing angle, for example, when the observer has to compare 8K UHD pictures. For those reasons, the number of pictures is limited to two for SDSCE or three for TSCES. This paper proposes a new method with multiple stimuli displayed in parallel on the same
. ' `
G # " % & % ' . ( % !- !) )!*+ & % % *- ! )!*+
!
monitor. A similar concept was introduced by ITU [3], where two full frame pictures were displayed on the same monitor, but this could work in the case of lower resolution video only. The proposed method does not have this restriction: the full resolution picture is segmented into an n-parts matrix and each cell can be Fig. 1. Diagram illustrating the method of splitting into an n-parts (n = 4) Rys. 1. Schemat ilustrujący metodę podziału ekranu na n-części (n = 4)
independently processed. The screen-splitting matrix defines the number and surface area of sub-images simultaneously given to the observer. The final image does not contain any lines separating the divisions (Fig. 1). The number n of stimuli can be changed according to experiments needs. The limit is the observers’ ability to perceive and evaluate more than one video at the same time.
2. Experiment To evaluate the proposed method, three experiments were conducted: PC, PC-1Ă—2, and PC-2Ă—2. The PC experiment was a pair comparison (PC) method in accordance with ITU recommendations [3]. The observers were to choose the video of a better quality in each pair of images presented alternately one after each other. In the PC-1Ă—2 experiment, each half of the frame was coded with different parameters and observers had to choose the better one. In the case of PC-2Ă—2, the screen was divided into two rows and two columns and observers were asked to point out which quarter had the best quality and which had the worst. The purpose of the PC experiment was to obtain reference results for PC-1Ă—2. The result of the PC-2Ă—2 experiment was supposed to answer the question of whether the observer was able to recognize the quality of an image divided into quarters, each of which was coded with different parameters. Thirty-nine observers participated in the PC experiment, 16 in the PC-1Ă—2 test, and 39 in the PC-2Ă—2 test. To allow a comparison of the obtained results, the same video was used in all experiments. This was a rendered animation composed of two characteristic scenes (Fig. 2). For both scenes, the
17
Subjective quality evaluation of a full resolution video with eligible number of reference pictures encoded with a certain bitrate was displayed at least once in any quarter of the screen
a)
v2 , In PC-2×2, the best part was always represented by a quarter encoded with 6 Mbps and the worst by one encoded with 3 Mbps. Therefore the difficulty of recognizing the quality of the image fragments in each set was constant. In the PC and PC-1×2 methods, the difficulty level varied for each set depending on the quality levels of the selected pair of images. Therefore, the results of the PC and PC-1×2 methods cannot be directly related to the results of PC-2×2. As expected, in the case of both the PC and PC-1×2 methods, the least votes were assigned to sequences encoded with the lowest bitrate (3 Mbps) and the most to the videos encoded with the highest bitrate (6 Mbps). The number of votes cast increases evenly with the increase of the bitrate, which means that in most cases observers were able to recognize and indicate an image of better quality. The increase in the number of votes is greater for PC-1×2, where two stimuli were presented simultaneously. The video with the lowest bitrate was chosen as the better one in only 2.6% of presentations in PC-1×2 and in only 4.9% of cases in the PC method. Assuming that choosing a video with a higher bitrate is the correct answer, the percentage of correct votes was calculated and is presented in Table 1 and Table 2. The correctness of PC-1×2 votes is higher than that achieved with the classical PC-1 method. Only in the pairs 3 & 5 Mbps and 5 & 6 Mbps there is a lower percentage of correct votes: − 3 & 5 Mbps pair: 93.75% (PC-1×2) and 97.44% (PC-1), − 5 & 6 Mbps pair: 81.25% (PC-1×2) and 82.05% (PC-1).
b)
Fig. 2. Frames of video test material: a) scene 1, b) scene 2 Rys. 2. Ramki materiału testowego: a) scena 1, b) scena 2
image was symmetrical to the centre of the screen to ensure the same probability of image distortion in each part. Due to the technical limitations of the laboratory equipment, a standard progressive MPEG-2 MP@HL (Main Profile @ High Level) with 1920 × 1080 px, 24 fps, and constant bitrate had to be chosen, but the concept is valid for any resolution and compression algorithm. The content was designed especially to obtain compression artefacts that are typical for MPEG-2: blocking, blurring, ringing, staircase, basis pattern artefacts, and mosquito noise [8]. Variation of the quality of the test sequences was achieved only by changing the bitrate in the coding process. Bitrate values were selected experimentally in pre-tests to provide clearly perceptible differences in quality. Four bitrate levels (n = 4) were chosen: 3, 4, 5, and 6 Mbps. Observers were not shown the source (uncompressed) video and were not informed about the compression parameters used. In the case of PC and PC-1×2, in accordance with recommendations [3], all n(n – 1) = 12 possible combinations were tested (Fig. 3). Because both sequences were displayed in parallel on the same screen in PC-1×2, the test time was halved in comparison with PC. In PC-2×2, the number of sequences evaluated simultaneously increased to four, therefore an observer would have to watch the same test sequence (with different sets of bitrates in individual quarters) 24 times. To lessen the observers’ fatigue and to shorten the duration of the experiment, 10 variations were chosen randomly (Fig. 4). In this test material, each sub-image
50 45 40 35 30 25 20 15 10 5 0
PC
345
PC-1x2
6 bitrate (Mbps)
Fig. 5. Distribution of votes cast on the better quality image in the PC and PC-1×2 methods Rys. 5. Rozkład głosów oddanych na sekwencję w eksperymentach PC i PC-1×2 w zależności od wielkości strumienia bitowego
5
3
6
4
4
5
5
6
6
5
3
6
5
4
6
3
4
3
3
4
4
6
3
5
Fig. 3. Stimulus presentation in the PC-1 and PC-1×2 study Rys. 3. Prezentacja materiału w eksperymencie PC-1 i PC-1×2
4 5
6 3
3 4
6 5
3 4
5 6
3 5
4 6
6 4
5 3
5 6
3 4
4 6
3 5
3 6
5 4
4 5
3 6
6 5
4 3
Fig. 4. Stimulus presentation in the PC-2×2 experiment Rys. 4. Prezentacja materiału w eksperymencie w eksperymencie PC-2×2
18
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
! "
5
6
3
x
69.23
97.44
94.87
4
87.17
x
92.31
100
5
94.87
92.31
x
82.05
6
97.44
94.87
74.36
x
Table 2. Percentage of correct votes in PC experiment PC-1×2 Tabela 2. Procentowy udział poprawnych odpowiedzi w eksperymencie PC-1×2
Mbps
3
4
5
6
3
x
81.25
93.75
100
4
93.75
x
100
100
5
100
100
x
81.25
6
100
100
87.5
x
In all other cases, the level of correctness of observers is higher for PC-1Ă—2. The number of pairs for which 100% of votes are correct has increased from one in the case of PC-1 to seven for PC-1Ă—2. The observers performed worst when the bitrate values were close to each other and close to the boundary values (highest or lowest) at the same time, namely the 3 & 4, 4 & 3, 5 & 6, and 6 & 5 Mbps pairs. The phenomenon was not observed for the 4 & 5 and 5 & 4 Mbps pairs. It is therefore likely that compression artefacts are still visible in the 4 Mbps video but disappear in the 5Â Mbps one. In contrast, visual differences in the cases of the 3 & 4, 4 & 3, 5 & 6, and 6 & 5 Mbps pairs were already harder to notice. In most cases, the PC-1Ă—2 test, where an observer viewed both images at the same time, gave a higher number of correct votes. This means that according to previous assumptions, the experiment was easier for participants and gave more accurate results. The results from the PC-2Ă—2 test are displayed in Fig. 6, which shows a general summary, for the whole study, of the total number of votes cast for specific bit streams considered to be the best (green) or the worst (red) quarters of the image. As the observer was watching all four bitrate representations at the same time, all votes cast for 6 Mbps as the best part and all votes cast for 3 Mbps as the worst may be indexed as correct. Therefore the percentage of votes cast correctly is very high: 80.51% in the case of the best quality and 87.95% in case of the
number of votes
400
87,95%
best
80,51%
worst
100%
100% correct answears
4
50% 0%
50%
top right
0% 36 bitrate (Mbps)
100%
36 bitrate (Mbps) 100% correct answears
3
correct answears
Mbps
worst. Next, in order of quality, 5 Mbps was denoted as the best by 18.46% of votes, while 4 Mbps was denoted as the worst by 11.28% of votes. A the same time there were two votes (0.51%) by observers indicating that the 4 Mbps and 3 Mbps streams were the best and three votes(0.77%) indicating that 5 Mbps streams were the worst. No observer voted for 6 Mbps as the worst. In general, it was easier for observers to identify quarters with the worst quality than those with the best. This shows, among other things, that humans focus primarily on visual distortions caused by excessive compression of the material. The better the image quality is, the more difficult it becomes to assess. The other research question was whether the correctness of the vote depends on the quarter in which the material is displayed. Figure 7 shows the number of correct indications for each quarter of the screen, taking into account indications of both the best and the worst parts. The correctness obtained in PC-2Ă—2 was high for each quarter. The lowest observed correctness level was 75.21%, in this case for the best quarter in the lower right quarter. The highest observed accuracy level was 92.31% for the worst part in the upper right quarter. In all quarters of the screen there was a consistent percentage of correct indications with an average value of 80.77% for the best quarter and 88.04% for the worst (average calculated for three quarters of the screen, as there was no 3 Mbps sequence in the lower left corner). None of the image quadrants were privileged by observers. The number of correct votes cast for the worst part was slightly higher than the number of correct votes cast for the best. This also shows once again that it is easier for observers to judge a decrease in image quality than an improvement. It may be assumed that there is no problem focusing attention
correct answears
Table 1. Percentage of correct votes in PC-1 experiment Tabela 1. Procentowy udział poprawnych odpowiedzi w eksperymencie PC-1
50% 0%
50% 0%
36 bitrate (Mbps)
36 bitrate (Mbps)
Fig. 7. Percentage of correct answers for the maximum (6 Mb/s) and the minimum (3 Mb/s) quality level depending on the frame quarter in which the stimulus was presented Rys. 7. Procentowy udział poprawnych odpowiedzi dla maksymalnego (6 Mb/s) oraz minimalnego (3 Mb/s) poziomu jakości w zaleşności od ćwiartki kadru, w którym prezentowany był bodziec
on the image divided into four parts and voting for two of them after viewing the test sequence
300 200 100
11,28% 0,51%
0,51%
x2 :
18,46% 0,77%
0%
0 3456 bitrate, Mbps Fig. 6. Percentage of votes given to a quarter with the best (green) and the worst (red) quality with respect to the bitrate Rys. 6. Procent głosów oddanych na ćwiartkę obrazu najlepszej (kolor zielony) i najgorszej (kolor czerwony) jakości w zaleşności od wielkości strumienia bitowego – większa wartość strumienia bitowego to lepsza jakość
Our idea of test material presentation made it possible to maintain more natural viewing conditions and to use a cheaper and more flexible experimental station. In the PC-1Ă—2 experiment, limiting the image quality information to half of the full screen did not disturb the image quality assessment, but the duration of the test was halved. Simultaneous presentation of the two stimuli resulted in increases in the number of scores cast correctly. What is more, the human audience was able to recognize and correctly evaluate four quarters coded with different bitrates displayed in parallel. Obtaining such promi-
19
Subjective quality evaluation of a full resolution video with eligible number of reference pictures sing results encourages further research involving experiments with natural videos instead of rendered ones and preferably of ultra-high resolution.
4. Kozamernik F., Steinmann V., Sunna P., Wyckens E., SAMVIQ—A New EBU Methodology for Video Quality Evaluations in Multimedia, “SMPTE Motion Imaging Journal�, Vol. 114, No. 4, 2005, 152–160, DOI: 10.5594/j11535. 5. ITU-R Recommendation BT.500-13: ‘Methodology for the subjective assessment of the quality of television pictures’, 2012 6. Hoffmann H., Itagaki T., Wood D., Hinz T., Wiegand T., A Novel Method for Subjective Picture Quality Assessment and Further Studies of HDTV Formats, “IEEE Transactions on Broadcasting�, Vol. 54, No. 1, 2008, 1–13, DOI: 10.1109/tbc.2008.916833. 7. Pinson M.H., Wolf S., Comparing subjective video quality testing methodologies, Visual Communications and Image Processing, 2003, DOI: 10.1117/12.509908. 8. Zeng K., Zhao T., Rehman A., Wang Z., Characterizing perceptual artifacts in compressed video streams, Proceedings of SPIE 9014, Human Vision and Electronic Imaging XIX, 90140Q (25 February 2014); DOI: 10.1117/12.2043128.
. # This work was fully supported by the statutory funds of Institute of Metrology and Biomedical Engineering, Warsaw University of Technology.
~ # ' 1. Shi Y., Sun H., Image and Video Compression for Multimedia Engineering, Image Processing Series, (CRC Press, Boca Raton, 2008), DOI: 10.1201/9781420007268. 2. Moldovan A.-N., Ghergulescu I., Muntean C. H., VQAMap: A Novel Mechanism for Mapping Objective Video Quality Metrics to Subjective MOS Scale, “IEEE Transactions on Broadcasting�, Vol. 62, No. 3, 2016, 610–627, DOI: 10.1109/tbc.2016.2570002. 3. ITU-T Recommendation P.910: Subjective video quality assessment methods for multimedia applications, 2008.
0 % % @ % F F < % Streszczenie: Artykuł prezentuje nową metodę oceny jakości materiału wideo. Koncepcja tej metody zakłada prezentację poddawanego ocenie materiału testowego równolegle z materiałem odniesienia, co ułatwia sformułowanie właściwej oceny. Nowatorskim pomysłem jest podział obrazu wideo w pełnej rozdzielczości na n-obrazową macierz. Kaşda komórka takiej macierzy moşe być zakodowana z innymi parametrami kompresji. W artykule zaprezentowano wyniki eksperymentów z zastosowaniem podziału na cztery części. Przeprowadzone analizy pokazują, şe moşliwe jest uzyskanie dokładniejszych wyników w krótszym czasie w porównaniu do klasycznej metody oceny. ' ` # & %
# # 3 !5
2 $ <= $ 3 !5
% '
' %
; % % O % # 8 % G F / # B & > # . ( 4 )!*- F % # % F # 6 ( % % % # % Z % ( F F # # F % %( % % # % % % 9 % = % ( H % Z = % % ( ?H % & % % % # 4. # 4 % % # % # % = # 4 F = % % = # % L % Z M& % 4. # = ( 9 % H & % ( %
; % % L)!!-M # > = # . ( L> .M = # 4 # B ( % A % ? ( ( # > . 9 = % [ % ( % % % Z
20
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
Pomiary Automatyka Robotyka, ISSN 1427-9126, R. 23, Nr 1/2019, 21–30, DOI: 10.14313/PAR_231/21
> # % %< ( % Robert Nebeluk, Piotr Marusak 0 > & > % ? . 4 # & X *E\*+& !!=JJE >
Streszczenie: Zastosowanie algorytmów regulacji predykcyjnej MPC do regulacji wielu procesów nieliniowych, o róşnym stopniu trudności, często umoşliwia osiągniecie bardzo dobrej jakości regulacji. Jest to moşliwe ze względu na odpowiednie uwzględnienie w strukturze algorytmu informacji uzyskanych z modelu procesu. Do formułowania zadania optymalizacji dla algorytmów regulacji predykcyjnej najczęściej zakładana jest stała trajektoria referencyjna dla całego horyzontu predykcji. W artykule przedstawiono moşliwości poprawy jakości regulacji przez zastosowanie trajektorii referencyjnej zmiennej na horyzoncie predykcji. Podczas porównywania jakości regulacji z wykorzystaniem trajektorii referencyjnych o róşnej postaci, wzięto pod uwagę dwie wielkości. Pierwszą z nich jest czas narastania wyjścia obiektu regulacji, a drugą – przeregulowanie. Badania prowadzono w układach regulacji dwóch obiektów: nieminimalnofazowego obiektu liniowego oraz nieliniowego reaktora polimeryzacji. Do regulacji w przypadku pierwszego obiektu, zastosowano algorytm DMC, a w przypadku drugiego – algorytm bazujący na modelu nieliniowym, z nieliniową predykcją i linearyzacją (NDMC-NPL). Przedstawione wyniki dobrze ilustrują moşliwości wpływania, za pomocą trajektorii referencyjnych o róşnych kształtach, na poprawę jakości regulacji oferowanej przez algorytmy predykcyjne. ( ` # & ( % & ( & F & I (
1. Wprowadzenie Algorytmy regulacji predykcyjnej MPC (ang. Model Predictive Control) są algorytmami, które wykorzystują model procesu do predykcji jego zachowania w przyszłości. W przypadku silnie nieliniowych obiektów regulacji, z duşym opóźnieniem, ograniczonymi sterowaniami lub wyjściami, zastosowanie algorytmów MPC zwykle przynosi dobrą jakość regulacji, co wynika z dostępu algorytmów do informacji o dynamice procesu. Algorytmy regulacji predykcyjnej równieş znają trajektorię zadaną wartości wyjść regulowanych oraz mają moşliwość wprowadzenia trajektorii referencyjnej, która jest w stanie modyfikować działanie algorytmu w celu poprawy jakości regulacji. Pojęcie trajektorii referencyjnej pojawiało się juş we wczesnych artykułach dotyczących regulacji predykcyjnej [10]. W niniejszym artykule ta ostatnia cecha algorytmów regulacji predykcyjnej będzie poddana szerszym badaniom. Algorytmy rozwaşane w artykule to wersje analityczne, mające mniejsze, w porów-
. ' `
0 B & 0 B ' % . ( % !- !) )!*+ & % % )K ! )!*+
!
naniu z algorytmami numerycznymi, moşliwości wpływania na zachowanie układu regulacji. Stąd wynika zainteresowanie autorów uşyciem mechanizmu trajektorii referencyjnej jako dodatkowej moşliwości strojenia takich algorytmów, rozszerzającej oferowane przez nie funkcjonalności. Znalezienie optymalnej trajektorii referencyjnej jest częstym problemem w robotyce. W artykule [4] jest to znalezienie trajektorii ruchu robota z punktu A do punktu B minimalizującej czas i wydatek energetyczny do wykonania przemieszczenia. Do wyznaczenia trajektorii referencyjnej jest uşyte odpowiednio sformułowane zadanie optymalizacji z ograniczeniami. Istnieją publikacje, w których generacja trajektorii referencyjnej dla robota jest wykonana na podstawie zdjęć otrzymanych przy pomocy kamery zainstalowanej na robocie [6] lub uporządkowanej listy punktów [9]. Rozwiązywane są równieş problemy związane z testowaniem utworzonego modelu robota przy wprowadzeniu wygenerowanej trajektorii referencyjnej [7]. Moşliwości wykorzystania wspomnianego mechanizmu w układach automatycznej regulacji liniowych obiektów są opisane w pracy [5]. Ponadto opisywana jest generacja trajektorii referencyjnej przez rozwiązanie problemu znalezienia najkrótszej ścieşki, stosując algorytm Dijkstry, którą później moşna zastosować w algorytmach MPC [8]. Podczas badań przetestowano róşne typy trajektorii referencyjnych. Eksperymenty przeprowadzono dla trajektorii o róşnym kształcie. W pierwszej kolejności zastosowano typowe trajektorie o postaci inercyjnej, a następnie – trajektorie nietypowe, skokowo zmienne, odpowiednio dobrane do kaşdego z obiektów.
21
? ) # M =#
Uşycie przedstawionych trajektorii ma za zadanie pokazanie moşliwości wpływania na cechy sygnału wyjściowego. Eksperymenty przeprowadzono dla obiektu liniowego oraz obiektu nieliniowego (reaktora polimeryzacji). W rozdziale 2 przedstawiono algorytmy regulacji predykcyjnej DMC oraz algorytm typu DMC bazujący na modelu nieliniowym, z nieliniową predykcją i linearyzacją (NDMC-NPL). W rozdziale 3 omówiono pierwszy obiekt regulacji – obiekt liniowy o trudnej dynamice oraz wyniki badań działania układu regulacji z zaimplementowanymi zmiennymi trajektoriami referencyjnymi. W rozdziale 4 przedstawiono nieliniowy obiekt regulacji (reaktor polimeryzacji) wraz z eksperymentami z wykorzystaniem trajektorii referencyjnych o róşnym kształcie, jak w rozdziale 3. Artykuł kończy podsumowanie otrzymanych wyników i wnioski z eksperymentów.
Przyjęto, şe czas narastania tn jest czasem potrzebnym, aby wyjście wzrosło od 0 do 100% swojej wartości końcowej [12].
Z2Y2 . # # ,-: Algorytm regulacji predykcyjnej DMC (ang. Dynamic Matrix Control) po raz pierwszy został zastosowany w przemyśle petrochemicznym i tam odniósł największy sukces. Regulator DMC uşywa liniowego modelu obiektu regulacji w postaci rzędnych odpowiedzi skokowej. Model ten opisany jest więc równaniem [1]:
(
)
yk +i |k =
i
∑s
n
n =1
⋅ Δuk −n +1 +
⎥y ⎤ ⎢ ⎼ ⎢y 2ref ⎼ =⎢ ⎼ ⎢ # ⎼ ⎢ ⎼ ⎢y ref ⎼ ⎣ N ⎌
⎥ s1 ⎢ ⎢ s2 M =⎢ ⎢# ⎢ ⎢s ⎣N
A
(6)
0
!
0
s1
!
0
#
%
#
⎤ ⎼ 0 ⎼ ⎼ # ⎼ ⎼ sN −Nu +1 ⎼⎌ 0
sN − 1 ! sN − N u + 2
(7)
OdpowiedĹş swobodna obiektu mozna obliczyc ze wzoru [1]: (8)
gdzie
(2)
R
(5)
y0 to odpowiedź swobodna, która zawiera wartości wyjścia obiektu otrzymane przy załoşeniu braku zmian sterowania na horyzoncie predykcji, y – przewidywane wartości wyjścia obiektu na horyzoncie predykcji.
⎥ s2 − s1 s3 − s2 ⎢ ⎢ s2 − s1 s 4 − s2 =⎢ ⎢ # # ⎢ ⎢s ⎣ N +1 − s1 sN + 2 − s2
! sD − 1 − sD − 2 sD − sD − 1 ⎤ ⎼ ! sD − sD − 2 sD − sD − 1 ⎼ ⎼ ⎼ % # # ⎼ ! sD − sD − 2 sD − sD −1 ⎼⎌
(9)
ΔuP – wektor przeszĹ‚ych przyrostĂłw sterowania: (3) ⎥ Δuk −1 ⎤ ⎢ ⎼ Δu P = ⎢ # ⎼ ⎢ ⎼ ⎢ Δu ⎼ ⎣ k − D +1 ⎌
gdzie: y(∞) – wartość wyjĹ›cia w stanie ustalonym, tm – czas wystÄ…pienia maksimum, czyli czas potrzebny, aby wyjĹ›cie osiÄ…gnęło pierwszy szczyt, y(tm) – maksymalna wartość wyjĹ›cia. I
⋅ Δuk −n +1 +sD ⋅ Δuk −D +i + dk
y 0 = y k + M P ⋅ Δu P
y(tm ) − y(∞) ⋅ 100% y(∞)
M
n
gdzie M × Δu – odpowiedź wymuszona obiektu regulacji zaleşna od przyszłych przyrostów sygnału sterującego Δu, M jest macierzą dynamiczną:
Kryteria porównawcze, według których będzie oceniana jakość regulacji, to przeregulowanie oraz czas narastania przebiegów wielkości regulowanych. Wielkość przeregulowania jest obliczana ze wzoru [11]:
O
∑s
n =i +1
(1)
MP
P
D −1
y = y 0 + M ⋅ Δu
N u −1 ⎪⎍ + Îť ∑ Δuk2+i |k ⎏ i =0 âŽâŽŞ
ref 1
22
(4)
gdzie dk = yk − ykM jest bĹ‚Ä™dem typu DMC i uwzglÄ™dnia niedokĹ‚adność modelowania oraz dziaĹ‚anie zakĹ‚ĂłceĹ„ niemierzalnych. Przewidywane wartoĹ›ci wyjĹ›cia moĹźna zapisać w postaci wektorowo-macierzowej [1]:
Dokonywana jest więc minimalizacja sumy przyszłych uchybów regulacji wyznaczonych na podstawie róşnicy między załoşoną trajektorią referencyjną a trajektorią przewidywanych wartości wyjść procesu, na horyzoncie predykcji, z uwzględnieniem kary za zmiany sygnału sterującego (drugi składnik minimalizowanego wskaźnika jakości). Podczas badań sprawdzano wpływ zmiennej trajektorii referencyjnej yref na horyzoncie predykcji na jakość regulacji. Trajektoria referencyjna będzie zatem wprowadzana do wskaźnika jakości jako następujący wektor:
yp =
⋅ Δuk −1 +sD ⋅ uk −D
Przewidywane wartości wyjścia obiektu są obliczane ze wzoru [1]:
gdzie Δu jest wektorem przyszłych (szukanych) przyrostów sterowania luk+i|k, yk+i|k jest przewidywaną w bieşącej chwili k wartością wyjścia obiektu dla przyszłej chwili k+i, wyznaczaną na podstawie modelu obiektu regulacji, l • 0 jest współczynnikiem waşącym przyszłe przyrosty sterowania, ykref+i |k to elementy trajektorii referencyjnej.
y ref
i
gdzie ykM – wyjście modelu regulacji w chwili k, si – rzędne odpowiedzi skokowej obiektu (i = 1, 2, ‌, D), D – horyzont dynamiki obiektu, który jest równy liczbie okresów próbkowania, po upływie których moşna odpowiedź obiektu uznać za ustaloną.
Algorytmy regulacji predykcyjnej w danej chwili k przewidują zachowanie procesu na wiele chwil do przodu, liczbę tych chwil nazywa się horyzontem predykcji N; natomiast ile zmian sterowania moşe nastąpić na tym horyzoncie określane jest jako horyzont sterowania Nu. Algorytmy MPC wykorzystują do predykcji model obiektu regulacji, stąd skróty nazw grup tych algorytmów. Przyszłe sterowania obliczane są na podstawie rozwiązywania zadania optymalizacji, czyli minimalizacji następującego wskaźnika jakości [1]: 2
∑s i =1
Z2 . # # '
⎪⎧ N min ⎨∑ ykref+i |k − yk +i |k Δu ⎊⎪ i =1
D −1
ykM =
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
(10)
A
NR 1 / 20 1 9
" $ % $ # #
yk – N-elementowy wektor wyjścia: ⎥yk ⎤ ⎢ ⎼ yk = ⎢ # ⎼ ⎢ ⎼ ⎢y ⎼ ⎣ k⎌
(11)
Analityczne rozwiązanie zadania optymalizacji (1) moşna przedstawić dla regulatora DMC następująco [1]:
(
Δu = M T ⋅ M + Ν ⋅ I
)
−1
(
⋅ M T ⋅ y ref − y 0
)
(12)
gdzie: I – macierz jednostkowa. W wyniku rozwiązania powyşszego zadania, otrzymany zostanie wektor przyszłych zmian sterowania Δu. Pierwszy element tego wektora słuşy do wyznaczenia sterowania uk, następnie optymalizacja jest powtarzana w kolejnej chwili próbkowania. Uwzględnienie ograniczeń w algorytmie analitycznym jest stosunkowo proste ze względu na moşliwość wykorzystania mechanizmu rzutowania sterowań na zbiór ograniczeń, który opiera się na następujących regułach [1]:
gdzie t jest opóźnieniem sygnału sterującego. Pierwszy element trajektorii swobodnej jest więc obliczany na podstawie aktualnego pomiaru wyjścia yk, wartości wyjścia i sterowania w przeszłych chwilach próbkowania oraz błędu typu DMC. Kolejne elementy są wyznaczane podobnie z tą róşnicą, şe zamiast wartości wyjścia w przyszłych chwilach próbkowania (z horyzontu predykcji), wykorzystywane są odpowiednie, juş wyznaczone elementy trajektorii swobodnej z poprzednich chwili z horyzontu predykcji. Poszczególne kroki algorytmu moşna przedstawić następująco: 1. Inicjalizacja; 2. Linearyzacja modelu nieliniowego; 3. Wyznaczenie macierzy dynamicznej na podstawie modelu zlinearyzowanego; 4. Obliczenie odpowiedzi swobodnej na podstawie modelu nieliniowego; 5. Sformułowanie zadania optymalizacji kwadratowej (1) oraz rozwiązanie go metodą analityczną, korzystając ze wzoru:
( )
Δu = ⎛⎜ M k �
T
−1
( ) â‹… (y
â‹… M k + Îť â‹… I ⎞âŽ&#x; â‹… M k âŽ
T
ref
− yk0
)
gdzie dla przyrostĂłw sterowania: − jeĹ›li uk|k < umin, to uk|k = umin; − jeĹ›li uk|k > umax, to uk|k = umax; oraz dla wartoĹ›ci sterowania: − jeĹ›li uk–1 + uk|k < umin, to uk|k = umin – uk–1; − jeĹ›li uk–1 + uk|k > umax, to uk|k = umax – uk–1.
⎥ yk0+1|k ⎤ ⎢ ⎼ yk0 = ⎢ # ⎼ ⎢ ⎼ ⎢y 0 ⎼ + k N | k ⎣ ⎌
6. Wyznaczenie nowego sterowania: W zaprezentowany sposób moşna łatwo zmodyfikować przyrosty sterowania generowane przez regulator predykcyjny.
uk = uk −1 + Δuk
(15)
Z2Z2 . # # $,-: $&9 W układzie regulacji nieliniowego reaktora polimeryzacji zastosowano, bazujący na modelu nieliniowym, algorytm typu DMC z nieliniową predykcją i linearyzacją (NDMC-NPL). Zasada działania algorytmu jest podobna do algorytmu DMC z tą róşnicą, şe w kaşdej iteracji dokonywana jest linearyzacja modelu nieliniowego, słuşąca do otrzymania macierzy dynamicznej, oraz wyznaczana jest odpowiedź swobodna obiektu na podstawie modelu nieliniowego. Macierz dynamiczna Mk zmienia się więc w kaşdej iteracji algorytmu i jest otrzymywana ze wzoru: ⎥ s1k ⎢ ⎢ s2k Mk = ⎢ ⎢# ⎢ ⎢s k ⎣N
0
!
0
s1k
!
0
#
%
#
sNk −1 ! sNk −Nu + 2
v2 ' ( #
% # 3.1. Opis obiektu Model obiektu liniowego opisany jest następującą transmitancją operatorową [2]:
G(s) =
⎤ ⎼ 0 ⎼ ⎼ # ⎼ ⎼ k sN −Nu +1 ⎼⎌ 0
(13)
Odpowiedź swobodna jest obliczana iteracyjnie w następujący sposób:
(
)
yk0+1|k = g yk , yk −1, ! , yk −nA , uk âˆ’Ď„ , uk −nB + dk y
0 k + 2|k
(
=g y
0 k + 1|k
)
, yk , ! , yk −nA +1 , uk âˆ’Ď„ +1 , uk −nB +1 + dk #
y
0 k + N |k
(
=g y
0 k + N −1|k
,y
0 k + N − 2|k
)
, ! , yk −nA + N −1, uk âˆ’Ď„ + N −1, uk −nB + N −1 + dk
(14)
(16)
Przez zastosowanie transformaty Z przy załoşonym okresie próbkowania Tp = 1 s, otrzymano model obiektu o postaci: G(z ) =
gdzie sik – rzędne odpowiedzi skokowej zlinearyzowanego w bieşącej chwili k obiektu.
−s + 1 e −4s 4s 2 + 2s + 1
−0,08427z + 0,277 −4 z z 2 − 1, 414z + 0,6065
(17)
Na rys. 1 przedstawiono pozyskaną odpowiedź skokową obiektu regulacji, którą wykorzystano w algorytmie DMC. Cechami charakterystycznymi obiektu regulacji jest opóźnienie t = 4 s, przejściowe dąşenie odpowiedzi w przeciwnym kierunku do kierunku dąşenia za wartością zadaną spowodowane dodatnim zerem oraz oscylacyjny charakter odpowiedzi. Ze względu na te właściwości obiekt jest stosunkowo trudny do regulacji mimo braku nieliniowości. W przypadku realizacji układu regulacji dla tego obiektu pomijany będzie błąd modelowania, zatem jako obiekt regulowany będzie przyjęty model liniowy, który wyznaczony został z dyskretnej postaci transmitancji modelu obiektu (17):
23
? ) # M =#
Czas narastania zwiększał się wraz ze wzrostem stałej czasowej T, ale malała równieş wielkość przeregulowania. Najmniejszą wielkość przeregulowania uzyskano dla największej stałej czasowej równej T = 1,6 s. Zauwaşmy, şe do pewnego stopnia jest moşliwe zmniejszenie wielkości przeregulowania przy braku zmiany czasu narastania, co moşna obserwować do wartości stałej czasowej T = 0,6 s. Następna rozwaşana trajektoria to inercja II rzędu:
(18)
Gref (s) =
K (T1s + 1)(T2s + 1)
(20)
W tej trajektorii K równieş jest równe wartości zadanej yzad oraz są dwie stałe czasowe T1 oraz T2, które będą zmieniane podczas badań. W tabeli 2 przedstawiono wyniki otrzymane dla tego przypadku. Badania dla tej trajektorii zaczęto od stałej czasowej T1 = 0,6 s ze względu na poprawę jakości regulacji, która miała miejsce przy omawianiu trajektorii (19). Zastosowanie T2 = 0,1 s zmniejszyło dodatkowo wielkość przeregulowania w porównaniu z zastosowaniem trajektorii (19). Czas narastania nie zmienił się. Rys. 1. Odpowiedź skokowa obiektu liniowego Fig. 1. Step response of the linear control plant
Tab. 2. Porównanie jakości regulacji dla róşnych parametrów trajektorii referencyjnej (20) Tab. 2. Comparison of control quality for different parameters of reference trajectory (20)
v2Z2 ~ '( (
%
*€ #
Rodzaj
Regulator DMC dla obiektu regulacji został zaimplementowany w wersji analitycznej, przedstawionej w rozdz. 2.1. Przyjęto następujące parametry regulatora: N = 20, Nu = 10, l = 2. Dobrane parametry regulatora umoşliwiły uzyskanie najkrótszego czasu ustalenia. Eksperymenty ze zmianą kształtu trajektorii referencyjnej będą przeprowadzone dla skoku wartości zadanej yzad = 10. Pierwsza rozwaşana postać trajektorii referencyjnej jest następująca:
Gref (s) =
K Ts + 1
(19)
Jest to inercja I rzędu, gdzie K jest równe róşnicy między wartością zadaną yzad a wartością początkową wyjścia obiektu y0, natomiast będzie prowadzona analiza zmiany jakości regulacji dla róşnych wartości stałej czasowej T. W tabeli 1 zamieszczono wyniki z tego badania.
trajektorii
yp [%]
tn [s]
Parametry
stała
3,3405
12
–
zmienna (20)
2,7161
12
T1 = 0,6 s, T2 = 0,1 s
zmienna (20)
2,5640
13
T1 = 0,6 s, T2 = 0,2 s
zmienna (20)
2,3696
13
T1 = 0,7 s, T2 = 0,1 s
zmienna (20)
2,9752
12
T1 = 0,5 s, T2 = 0,1 s
zmienna (20)
2,8588
12
T1 = 0,5 s, T2 = 0,2 s
zmienna (20)
2,6831
13
T1 = 0,5 s, T2 = 0,3 s
Następnie przeprowadzono eksperymenty dla pierwszej nietypowej postaci trajektorii, w celu zmniejszenia czasu narastania:
Tab. 1. Porównanie jakości regulacji dla róşnych parametrów trajektorii referencyjnej (19) Tab. 1. Comparison of control quality for different parameters of reference trajectory (19) Rodzaj trajektorii
yp [%]
tn [s]
Parametry
stała
3,3405
12
–
zmienna (19)
3,3404
12
T = 0,1 s
zmienna (19)
3,0482
12
T = 0,5 s
zmienna (19)
2,8201
12
T = 0,6 s
zmienna (19)
2,5083
13
T = 0,7 s
zmienna (19)
1,81
13
T=1s
zmienna (19)
0,0679
24
T = 1,6 s
24
P
O
M
I
A
R
Y
•
A
U
T
O
⎧⎪K ⋅ y zad dla k < m y(ref k) = ⎨ ⎪⎊ y zad dla k ≼ m
(21)
W tym przypadku K jest parametrem nałoşonym na wartość zadaną dla określonych chwil k. W tej trajektorii będą zmieniane parametry K oraz m, gdzie m jest zmienną określającą zakres elementów, które będą zmieniane w zaleşności od wartości parametru K. Wyniki uzyskane przy zastosowaniu tej postaci trajektorii referencyjnej zostały zamieszczone w tabeli 3. Celem tych testów było sprawdzenie, jak zmiana wartości oraz liczby zmienionych pierwszych elementów trajektorii referencyjnej wpłynie na jakość regulacji. Dla m = 2 i K = 0,1 uzyskano wydłuşenie czasu narastania przy jednoczesnym zmniejszeniu wielkości przeregulowania. W przypadku zwiększenia wartości parametru K uzyskano lepszą jakość regulacji, więc kolejne eksperymenty wykonano ze zmienianym parametrem m. Dla m = 3 przeregulowane było najmniejsze, przy niezmienionym czasie narastania, więc przeprowadzono dla tego parametru dodatkowy eksperyment dla K > 1. M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
" $ % $ # #
Tab. 3. Porównanie jakości regulacji dla róşnych parametrów trajektorii referencyjnej (21) Tab. 3. Comparison of control quality for different parameters of reference trajectory (21) Rodzaj trajektorii
yp [%]
tn [s]
Parametry
stała
3,3405
12
zmienna (21)
3,2324
zmienna (21)
Tab. 4. Porównanie jakości regulacji dla róşnych parametrów trajektorii referencyjnej (22) Tab. 4. Comparison of control quality for different parameters of reference trajectory (22) Rodzaj trajektorii
yp [%]
tn [s]
Parametry
–
stała
3,3405
12
–
13
K = 0,1, m = 2
zmienna (22)
0,4541
16
K1 = 0,9, K2 = 0,1, m = 3
3,2245
12
K = 0,9, m = 2
zmienna (22)
1,1660
13
K1 = 0,9, K2 = 0,8, m = 3
zmienna (21)
2,6817
12
K = 0,9, m = 3
zmienna (22)
2,8476
12
K1 = 1,1, K2 = 0,9, m = 3
zmienna (21)
1,7089
13
K = 0,9, m = 4
zmienna (22)
3,7466
11
K1 = 1,2, K2 = 0,9, m = 3
zmienna (21)
3,9993
12
K = 1,1, m = 3
zmienna (22)
2,3548
12
K1 = 1,2, K2 = 0,8, m = 3
zmienna (22)
1,8989
11
K1 = 1,8, K2 = 0,2, m = 3
zmienna (22)
2,5422
10
K1 = 1,9, K2 = 0,1, m = 3
zmienna (22)
0,1604
20
K1 = 1,8, K2 = 0,1, m = 3
Dla tej trajektorii, mimo moşliwości dodatkowego zmniejszenia wielkości przeregulowania dla parametrów K = 0,9 i m = 3 w porównaniu z zastosowaniem trajektorii (20), nie udało się skrócić czasu narastania. Nietypowa postać trajektorii pokazuje jednak dalsze moşliwości poprawy jakości regulacji. Ostatnia badana trajektoria referencyjna to:
y(ref k)
⎧K 1 ⋅ y zad dla k < m ⎪ ⎪ = ⎨K 2 ⋅ y zad dla k = m ⎪ zad dla k > m ⎪⎊ y
(22)
Podczas badań nad trajektorią (21) zmieniano wielkość przewidywanych uchybów regulacji z pierwszych chwil działania algorytmu. Obecnie rozwaşana jest postać trajektorii, w której następuje zwiększenie przewidywanych uchybów regulacji w określonym przedziale czasu, a następnie zmniejszenie uchybu regulacji w pojedynczej chwili z horyzontu predykcji. Zmieniane parametry trajektorii to K1 i K2. Przyjęto parametr m = 3, w związku z poprawą jakości regulacji przy eksperymentach z trajektorią (21). W tabeli 4 przedstawiono wyniki badań po zastosowaniu trajektorii (22). Podczas badań z wykorzystaniem tej trajektorii starano się wpłynąć korzystnie na czas narastania przez nadanie większej wartości parametrowi K1, w porównaniu z wartością parametru K2. W związku z tym, şe dla trajektorii (21) uzyskano lepsze rezultaty przy zwiększaniu wartości parametru K, to analogicznie postanowiono zwiększać parametr K1 trajektorii (22). Dla K1 = 1,2 otrzymano pierwszy wynik ze skróceniem czasu narastania. Okazało się, şe dopiero odpowiednio duşe K1 i odpowiednio małe K2 daje najlepszą jakość regulacji, równieş pod względem redukcji przeregulowania. Niestety skrócenie czasu narastania moşe odbywać się kosztem zwiększenia przeregulowania, dlatego dobór parametrów powinien być wykonywany stopniowo w kolejnych eksperymentach. Najkrótszy czas narastania, równy tn = 10 s, uzyskano dla K1 = 1,9, K2 = 0,1. Wielkość przeregulowania przy zastosowaniu tej trajektorii równieş była mniejsza niş przy zastosowaniu stałej trajektorii referencyjnej. Najmniejsza wielkość przeregulowania została otrzymana dla parametrów K1 = 1,8, K2 = 0,1. Zastosowanie typowych i nietypowych kształtów trajektorii referencyjnych dla obiektu liniowego pokazało moşliwości wpływania na czas narastania oraz przeregulowanie. Najskuteczniejszymi trajektoriami pod względem zmniejszenia wielkości przeregulowania są typowe trajektorie inercyjne (19) i (20).
Natomiast jest to rĂłwnieĹź moĹźliwe przy zastosowaniu trajektorii nietypowych jak (22). Ponadto trajektorie typu (22) umoĹźliwiajÄ… skutecznÄ… redukcjÄ™ czasu narastania.
x2 ' ( #
' 4.1. Opis obiektu Drugim obiektem regulacji jest reaktor chemiczny, w którym zachodzi proces polimeryzacji, czyli zamiany związków chemicznych o małej masie cząsteczkowej (monomery) na związki o większej masie cząsteczkowej (polimery). Uzyskiwane polimery są powszechnie stosowane w produkcji tworzyw sztucznych, a takşe produktów chemicznych takich jak: farby, lakiery, oleje przemysłowe, środki smarujące czy kleje. Reaktor opisany jest następującymi równaniami [3]:
x 1 = 10 6 − x1 − 2,4568x1 x 2
(
)
(23)
x 2 = 80u − 10,1022x 2
(24)
x 3 = 0,0024121x1 x 2 + 0,112191x 2 − 10x 3
(25)
x 4 = 245,978x1 x 2 − 10x 4
(26)
a równanie wyjścia obiektu jest dane jako:
y=
x4 x3
(27)
Wartości poszczególnych zmiennych w punkcie pracy są następujące [3]:
25
? ) # M =#
x10 = 5,50677 kmol/m3, x20 = 0,132906 kmol/m3,
Dla przedstawionego nieliniowego obiektu regulacji rĂłwnieĹź przeprowadzono eksperymenty z zastosowaniem trajektorii typowych (inercyjnych) jak i nietypowych. OprĂłcz trajektorii (19), (20) i (22), zastosowano rĂłwnieĹź trajektoriÄ™ o postaci:
x30 = 0,0019752 kmol/m3, x40 = 49,3818 kmol/m3, u0 = 0,016783 m3/h, y0 = 25000,5.
⎧⎪K ⋅ y zad dla k ≤ m = y(ref ⎨ zad k) dla k > m ⎪⎊ y
Zadanie regulacji polega na sterowaniu przepływem inicjującym monomerów u tak, aby regulować cięşar molowy produktu y. Charakterystyka statyczna obiektu regulacji jest przedstawiona na rys. 2 i ilustruje nieliniowy charakter obiektu.
(28)
Parametrami trajektorii są K i m, które pełnią tę samą rolę w projektowaniu trajektorii, jaką pełnią w przypadku trajektorii (21). Eksperymenty przeprowadzono teş dla trajektorii o postaci wyznaczonej na podstawie trajektorii (22):
y(ref k)
x2Z2 ~ '( (
%
*€ #
Dla trajektorii (19) zmniejszenie przeregulowania wiąşe się ze wzrostem czasu narastania, tak jak w przypadku poprzedniego obiektu regulacji. Zaczynając od stałej czasowej T = 0,02 i zwiększając ją, zmniejszano wartość przeregulowania y p1 i y p2 wraz ze zwiększaniem czasów narastania tn1 i tn2 . Dla T = 0,5 uzyskano brak przeregulowania. Przebiegi wyjścia i sterowania z zastosowaniem trajektorii (19) z parametrem T = 0,5 przedstawiono na rys. 4. Podobnie jak w przypadku poprzedniego obiektu regulacji, trajektorie inercyjne najskuteczniej minimalizowały wielkość przeregulowania. Dla trajektorii (20), jak w badaniach nad obiektem liniowym, przyjęto T1 = 0,02, gdzie czas narastania był prawie taki sam, jak przy zastosowaniu stałej trajektorii, a przeregulowanie było mniejsze. Przystąpiono do stopniowego zwiększania drugiej stałej czasowej T2. Otrzymano wydłuşenie czasu narastania oraz niewielkie zmniejszenie przeregulowania. W przeciwieństwie do obiektu liniowego nie ma wyraźnej poprawy jakości regulacji w wyniku zastosowania tej trajektorii. W przypadku trajektorii (28) istotny był dobór parametrów m i K. Na podstawie badań nad obiektem liniowym, w początkowych eksperymentach przyjęto wartość parametru K bliską 1. Parametr m zwiększano stopniowo i analizowano zmieniającą się jakość regulacji. W związku z tym, şe okres próbkowania Tp = 0,01, to dopiero dla m = {25, 50} było widoczna róşnica w jakości regulacji. Wybrano m = 50 do dalszych badań, bo powodowało zadowalające zmniejszenie wielkości przeregulowania kosztem zwiększenia czasu narastania. Dla m = 100 uzyskano największą redukcję wielkości przeregulowania. Stosując trajektorię (22), podobnie jak przy poprzednim obiekcie regulacji, próbowano skompensować zwiększenie przewidywanych uchybów regulacji przez parametr K1 przez odpowiednie zmniejszenie uchybu z chwili określonej przez parametr m, wykorzystując parametr K2 tak, aby uzyskać krótszy czas narastania. Dla mniejszych wartości parametru K2 czas narastania wydłuşał się. Najkrótszy czas narastania otrzymano dla tej trajektorii przy m = 50, K1 = 1,1 i K2 = 0,9. Stopniowe zwięk-
Dla przedstawionego obiektu regulacji zaprojektowano regulator NDMC-NPL. Dobrano następujące wartości parametrów: N = 50, Nu = 10, l = 2e12, oferujące najkrótszy czas ustalenia przebiegu wyjściowego. Załoşono okres próbkowania Tp = 0,01 h. Dla tak dobranych parametrów dokonano porównania działania regulatorów. Ze względu na uzyskany najkrótszy czas ustalenia oraz najmniejszą wielkość przeregulowania, do eksperymentów z róşnymi kształtami trajektorii referencyjnych, wybrano regulator NDMC-NPL zaimplementowany w wersji analitycznej z uwzględnionym ograniczeniem dolnym wartości sterowania umin = 0. Na rys. 3 przedstawiono odpowiedź skokową uzyskaną po linearyzacji modelu opisanego równaniami (23)–(27). Przedstawiona odpowiedź skokowa obiektu ma charakter inercyjny i ma stosunkowo krótki czas ustalenia.
Rys. 3. Odpowiedz skokowa zlinearyzowanego modelu obiektu Fig. 3. Step response of the linearized polymerization reactor model
P
O
M
I
A
R
Y
•
A
U
(29)
W tym przypadku parametry trajektorii K1 i K2 zmieniają wielkość przewidywanego uchybu w zaleşności od parametru m. Parametr m w przeciwieństwie do trajektorii (22) definiuje przedział czasowy, w którym dwa przewidywane uchyby regulacji mają mniejszą wartość od wcześniejszych. Wyniki badań dla opisanych trajektorii zostały przedstawione w tab. 5. Podano w niej wartość czasu narastania i przeregulowania dla wartości zadanych najbardziej oddalonych od punktu pracy. Wielkości y p1 , tn1 zatem dotyczą przebiegów dla wartości zadanej równej 29 500,5 a y p2 , tn2 – dla wartości zadanej równej 20 500,5.
Rys. 2. Charakterystyka statyczna reaktora polimeryzacji Fig. 2. Steady-state characteristic of the polymerization reactor
26
⎧ K 1 ⋅ y zad dla k < m ⎪ ⎪ = ⎨K 2 ⋅ y zad dla k = {m, m + 1} ⎪ y zad dla k > m + 1 ⎪⎊
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
" $ % $ # #
Tab. 5. Porównanie jakości regulacji przy zastosowaniu trajektorii referencyjnej o róşnym kształcie Tab. 5. Comparison of control quality for different shapes of reference trajectories
Rodzaj trajektorii
y p1 [%]
y p2 [%]
tn1 [h]
tn2 [h]
Parametry
stała
0,4203
0,1141
1,31
1,92
–
(19)
0,4201
0,1122
1,32
1,94
T = 0,02
(19)
0,4192
0,1006
1,35
2,01
T = 0,05
(19)
0,4131
0,0597
1,41
2,23
T = 0,10
(19)
0,1928
0
1,72
–
T = 0,30
(19)
0
0
–
–
T = 0,50
(20)
0,4201
0,1120
1,33
1,95
T1 = 0,02, T2 = 0,005
(20)
0,4201
0,1116
1,33
1,95
T1 = 0,02, T2 = 0,01
(20)
0,4200
0,1100
1,34
1,97
T1 = 0,02, T2 = 0,02
(20)
0,4190
0,0985
1,37
2,04
T1 = 0,02, T2 = 0,05
(28)
0,4202
0,1141
1,32
1,93
K = 0,1, m = 2
(28)
0,4203
0,1141
1,31
1,92
K = 0,9, m = 2
(28)
0,4203
0,1140
1,31
1,92
K = 0,9, m = 3
(28)
0,4203
0,1137
1,31
1,93
K = 0,9, m = 4
(28)
0,4203
0,1134
1,31
1,93
K = 0,9, m = 5
(28)
0,4199
0,1103
1,32
1,94
K = 0,9, m = 10
(28)
0,4174
0,0835
1,33
2,05
K = 0,9, m = 25
(28)
0,3978
0,0046
1,37
2,99
K = 0,9, m = 50
(28)
0,1787
0
1,62
–
K = 0,9, m = 100
(22)
0,4596
0,3667
1,24
1,51
K1 = 1,1, K2 = 0,9, m = 50
(22)
0,4238
0,1521
1,28
1,81
K1 = 1,1, K2 = 0,9, m = 25
(22)
0,4341
0,2642
1,26
1,61
K1 = 1,1, K2 = 0,1, m = 50
(29)
0,4076
0,1627
1,28
1,76
K1 = 1,1, K2 = 0,1, m = 50
(29)
0,4566
0,5197
1,21
1,37
K1 = 1,2, K2 = 0,1, m = 50
(29)
0,4049
0,1525
1,29
1,78
K1 = 1,1, K2 = 0,05, m = 50
(29)
0,4526
0,5034
1,21
1,38
K1 = 1,2, K2 = 0,05, m = 50
27
? ) # M =#
Rys. 4. Odpowiedzi układu regulacji z regulatorem NDMC-NPL na zmiany wartości zadanej y zad = {20 500,5; 22 000,5; 23 500,5; 25 000,5; 26 500,5; 28 000,5; 29 500,5}; z zastosowaniem trajektorii referencyjnej (19) dla T = 0,5 – linia ciągła; z zastosowaniem stałej trajektorii referencyjnej – linia przerywana; góra – wyjście, dół – sterowanie Fig. 4. Responses of the control system with NDMC-NPL algorithm to setpoint changes to yzad = {20 500.5; 22 000.5; 23 500.5; 25 000.5; 26 500.5; 28 000.5; 29 500.5}; with reference trajectory (19) for T = 0.5 – solid line, with constant reference trajectory – dashed line; above – output, below – control
Rys. 5. Odpowiedzi układu regulacji z regulatorem NDMC-NPL na zmiany wartości zadanej y zad = {20 500,5; 22 000,5; 23 500,5; 25 000,5; 26 500,5; 28 000,5; 29 500,5}; z zastosowaniem trajektorii referencyjnej (29) dla K1 = 1,2 i K 2 = 0,1 – linia ciągła; z zastosowaniem stałej trajektorii referencyjnej – linia przerywana; góra – wyjście, dół – sterowanie Fig. 5. Responses of the control system with NDMC-NPL algorithm to set– point changes to yzad = {20 500.5; 22 000.5; 23 500.5; 25 000.5; 26 500.5; 28 000.5; 29 500.5}; with reference trajectory (29) for K1 = 1.2 and K2 = 0.1 – solid line, with constant reference trajectory – dashed line; above – output, below – control
szanie parametru K1 i jednoczesne zmniejszanie parametru K2 doprowadziło więc do osiągniecia celu. Celem implementacji trajektorii (29) było uzyskanie dodatkowego zmniejszenia czasu narastania przez zmianę liczby przewidywanych uchybów regulacji objętych wpływem zmiany parametru K2. W sumie dwa uchyby były objęte tą zmianą. Podobnie jak w trajektorii (22) stopniowo zwiększano parametr K1 i jednocześnie zmniejszano parametr K2. Wykonane eksperymenty doprowadziły do otrzymania wartości parametrów, dla których został uzyskany najkrótszy czas narastania. Dla K1 = 1,2 i K2 = 0,1 otrzymano najkrótszy czas narastania, jednak z większym przeregulowaniem w porównaniu ze stałą trajektorią referencyjną. Przebiegi wyjścia i sterowania z zastosowaniem trajektorii (29) oraz parametrów K1 = 1,2 i K2 = 0,1, przedstawiono na rys. 5. Dla nieliniowego obiektu regulacji zastosowanie typowych (inercyjnych) i nietypowych trajektorii referencyjnych wskazało wiele moşliwości wpływania na jakość regulacji. W przypadku minimalizacji wielkości przeregulowania, trajektorie
inercyjne dawały najlepsze rezultaty. Dzięki trajektoriom o nietypowym kształcie moşliwe było natomiast skuteczne skracanie czasu narastania. Aby uzyskać jeszcze lepszą jakość regulacji, zwłaszcza pod kątem skrócenia czasu narastania, moşna stosować nietypowe trajektorie referencyjne i dalej modyfikować trajektorię (29) do bardziej skomplikowanych postaci. Moşna m.in. rozwaşyć zwiększanie zakresów wyznaczanych przez parametr m lub dodawanie nowych.
28
P
O
M
I
A
R
Y
•
A
U
T
O
M
X2 & Podsumowując, w artykule zostały przeprowadzone badania z zastosowaniem róşnych kształtów trajektorii referencyjnych w układach regulacji predykcyjnej dwóch obiektów: nieminimalnofazowego obiektu liniowego oraz nieliniowego reaktora polimeryzacji. Dla obu obiektów przetestowano trajektorie inercyjne, które wspominane są w literaturze [1, 5] jako trajektorie mogące wpływać na działanie układów regulacji. Ponadto A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
" $ % $ # #
zostały zaproponowane i przetestowane trajektorie o nietypowym kształcie. Badania wykazały, şe trajektorie inercyjne są najprostsze do wyznaczenia i stosując je moşna skutecznie zredukować wielkość przeregulowania, kosztem wydłuşenia lub braku zmiany czasu narastania, w zaleşności od przyjętych stałych czasowych. Aby poprawić jakość regulacji naleşy stopniowo zwiększać jedną ze stałych czasowych i analizować wpływ zmiany na obie wielkości. Przeregulowanie będzie wraz ze wzrostem stałej czasowej maleć. W chwili wydłuşenia się czasu narastania naleşy zaprzestać tego działania. W ten sposób zostanie wyznaczona wartość stałej czasowej, dla której otrzymamy mniejszą wartość przeregulowania przy braku zmiany czasu narastania. Zastosowanie trajektorii nietypowych, jak opisane w artykule lub podobnych, moşe przynieść skrócenie czasu narastania lub wielkości przeregulowania. Dla trajektorii (21) naleşy przyjąć parametr K bliski wartości 1, aby zmniejszyć wielkość przeregulowania bez zmiany czasu narastania. W przypadku trajektorii (22), (28) i (29) naleşy stopniowo zwiększać parametr K1 oraz jednocześnie zmniejszać parametr K2, aby znaleźć wartości parametrów, dla których zostanie otrzymany krótszy czas narastania.
Bibliografia 1. Tatjewski P., Sterowanie zaawansowane obiektów przemysłowych. Struktury i algorytmy. Wydanie drugie zmienione., Akademicka Oficyna Wydawnicza EXIT, Warszawa 2016. 2. Marusak P., Regulatory predykcyjne z załoşoną trajektorią przyrostów sterowania i uwzględnianiem ograniczeń sygnału sterującego, „Pomiary Automatyka Robotyka�, R. 12, Nr 2, 2008, 581–590. 3. Doyle F.J., Ogunnaike B.A., Pearson R.K., Nonlinear Model-based Control Using Second order Volterra Models., „Automatica�, Vol. 31, No. 5, 1995, 697–714.
4. Ardakani M.G., Olofsson B., Robertsson A., Johansson R., Real-Time Trajectory Generation using Model Predictive Control, „2015 IEEE Conference on Automation Science and Engineering(CASE)�, 2015, 942–948. 5. Seborg D.E., Edgar T.F., Mellichamp D.A., Process dynamics and control, John Wiley & Sons, 2011. 6. Burlacu A., Lazar C., Reference Trajectory-Based Visual Predictive Control, „Advanced Robotics�, Vol. 26, No. 8–9, 2012, 1035–1054. 7. Krasnanský R., Valach P., Soós D., Zarbakhsh J., Reference trajectory tracking for a multi-DOF robot arm, „Archives of Control Sciences�, Vol. 25, No. 4, 2015, 513– 527, DOI: 10.1515/acsc-2015-0033. 8. Wu Q., Xi Y., Nagy Z., Li D., Economic optimization in transient processes for model predictive control with a dynamic reference trajectory, „Computers and Chemical Engineering�, Vol. 121, 2019, 224–231, DOI: 10.1016/j.compchemeng.2018.11.006. 9. Haghighat S., Di Cairano S., Bortoff S.A., Method for computing reference trajectory to position machine, 2015. 10. Richalet J., Rault A., Testud J.L., Papon J., Model algorithmic control of industrial processes, „IFAC Proceedings Volumes�, Vol. 10, No. 16, 1977, 103–120, DOI: 10.1016/S1474-6670(17)69513-2. 11. Ahmetzyanov I.Z., Dem’yanov D.N., Determination of the Overshoot Scalar Control Systems with Transfer Zero and Binomial Law of Poles Distribution, „Indian Journal of Science and Technology�, Vol. 10, No. 1, 2017, DOI: 10.17485/ijst/2017/v10i1/109928. 12. Aly A.A., Salem F.A., A New Accurate Analytical Expression for Rise Time Intended for Mechatronics Systems Performance Evaluation and Validation, „International Journal of Automation, Control and Intelligent Systems�, Vol. 1, No. 2, 2015, 51–60.
4 ] ( 0 % 0 # F G # . ( Abstract: Applying model predictive control (MPC) algorithms to control many processes, of different difficulty level, often allows improving control quality. It is possible by including information received from a process model in the algorithm structure. When defining the optimization problem for the predictive control algorithms most often a time-constant reference trajectory is assumed. Possibilities of improving the control quality by applying a time reference trajectory variable on the prediction horizon are presented in the paper. Two quantities are considered when comparing control quality: the rise time of the control plant output, and the overshoot. The experiments were conducted in the control systems of two control plants: a linear nonminimumphase plant and a nonlinear polymerization reactor. In the control system of the first control plant the DMC predictive control algorithm was used. For the nonlinear reactor the NDMC-NPL algorithm based on a nonlinear model was applied. It is demonstrated that by using the reference trajectories of different shapes it is possible to improve control quality offered by the MPC control algorithms. Keywords` # & % & & & Z
29
? ) # M =#
"
"
G X F ' % %
0 B ' % % > *+Q- 6 ( = @ *++Q % )!! & F % > % ? . 4 # 0 = > 0 4 8 = 4 # 0 > R ( F S[ F < & ( ( % & % & % ( % % & % = S
% > *++) 6 @ )!*E > % 4 @ 0 % 0 > = & ( @ N )!*+ > % ? . 4 # = 0 > & F % 8 G F 0 % 8 4 @ $ = ( 4 8 4 # 0 = > R ( [ ( = ( < & (< ( ( = % ( % % (
30
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
Pomiary Automatyka Robotyka, ISSN 1427-9126, R. 23, Nr 1/2019, 31–42, DOI: 10.14313/PAR_231/31
> F ^ < ( ( S ( % ( 046 ) ;" $ >
3 4 $
3 $ # 5" 5 >
0 > & > % ? . 4 # & 4 8 4 # & X *E\*+& !!=JJE >
4 Praca ma na celu zbadanie i porównanie regulatorów PID o parametrach uzyskanych w wyniku optymalizacji wybranych wskaźników jakości pod kątem właściwości generowanych sygnałów sterujących. Punktem wyjścia do analizy są symulacje przeprowadzone w środowisku MATLAB przeprowadzone dla pięciu obiektów (z czterech klas) na sześciu typach wskaźników. Przedstawiono w szczegółach zastosowane metody i wykorzystane algorytmy. W pracy prezentowane są otrzymane w trakcie optymalizacji nastawy regulatorów PID, przebiegi sygnałów procesowych w badanych układach regulacji oraz obliczone parametry sygnałów sterujących, a takşe sformułowane na podstawie badań obserwacje i wnioski. ( ` ( 046& & & ^ & ( S
1. Wprowadzenie W związku z postępującym rozwojem techniki i rosnącą w wyniku globalizacji konkurencją coraz większą uwagę zwraca się na dokładne sterowanie procesami ciągłymi, które stanowią podstawę wielu gałęzi przemysłu. Przejawia się to rosnącą popularnością zaawansowanych algorytmów regulacji, takich jak sterowanie predykcyjne [22], jednak jego stosowalność ma swoje praktyczne ograniczenia [19] i najczęściej stosowanym rozwiązaniem pozostaje występujący w róşnych formach regulator PID [12]. Regulacja z wykorzystaniem tego algorytmu ma wiele zalet, które przyczyniły się do jego popularności w przemyśle. Regulatory tego rodzaju nie wymagają modelu obiektu do działania, nie wymagają duşego nakładu obliczeń ze strony sterowników, a często mimo swojej prostoty zapewniają w zupełności wystarczającą dokładność i szybkość regulacji [11]. W celu zapewnienia poprawnego działania układu regulacji, nastawy (parametry) regulatora PID powinny zostać odpowiednio dobrane pod kątem właściwości obiektu sterowania. Poprawne dostrojenie poprawia jakość produktu, ogranicza koszty eksploatacji, a takşe podnosi bezpieczeństwo obiektu.
. ' `
0 6 6 T & % ' % . ( % - !) )!*+ & % % )E ! )!*+
!
Trafny dobĂłr nastaw regulatora ma istotny wpĹ‚yw na wiele aspektĂłw procesu, ale przede wszystkim prowadzi do zwiÄ™kszenia zyskĂłw i ograniczenia kosztĂłw, co powinno być wystarczajÄ…cym argumentem za tym, Ĺźe jest to istotny problem, ktĂłry nie powinien być lekcewaĹźony [5, 13]. W literaturze temat ten rozpatrywany jest od wielu lat. Pierwsze raportowane prace pojawiajÄ… siÄ™ w latach 60. wraz z ocenami pracy ukĹ‚adĂłw sterowania w papiernictwie dokonanymi przez Ă…strĂśma [1]. Badania byĹ‚y kontynuowane w kolejnych dziesiÄ™cioleciach [23, 26]. Znaczne przyspieszenie nastÄ…piĹ‚o w 1989 r. wraz z pracami Harrisa [9] na temat metod minimalno-wariacyjnych. Obecne rozwiÄ…zania pozwalajÄ… na zastosowanie dla róşnych przypadkĂłw oraz zjawisk, jak na przykĹ‚ad nieliniowoĹ›ci [10], ukĹ‚adĂłw wielkiej skali [14], regulacji predykcyjnej [16], systemĂłw wielowymiarowych [27], uwzglÄ™dniajÄ…c wpĹ‚yw tarcia w urzÄ…dzeniach wykonawczych [4], zjawiska oscylacji [20] i wiele innych [8]. Owe badania znalazĹ‚y rĂłwnieĹź przeĹ‚oĹźenie komercyjne, jako Ĺźe istnieje wiele rozwiÄ…zaĹ„ informatycznych automatyzujÄ…cych proces [3]. PoniewaĹź opracowano wiele róşnych wskaĹşnikĂłw, moĹźna wydzielić kilka grup podejść [8, 11, 21]: − rozwiÄ…zania w dziedzinie czasu wykorzystujÄ…ce odpowiedĹş skokowÄ…: uchyb ustalony, przeregulowanie, czas regulacji, dojĹ›cia, osiÄ…gniÄ™cia maksimum, wskaĹşniki: Area Index, Output Index, R-index, Idle Index, − wskaĹşniki bazujÄ…ce na przebiegach czasowych zmiennych biorÄ…cych udziaĹ‚ w sterowaniu (sygnaĹ‚ sterujÄ…cy, uchyb regulacji, wyjĹ›cie procesu), najczęściej w postaciach caĹ‚kowych: Ĺ›redni bĹ‚Ä…d kwadratowy MSE (Mean Square Error), Ĺ›rednia caĹ‚ka z moduĹ‚u IAE (Integral of Absolute Error), ale i innych jak wskaĹşnik amplitudowy AMP (Amplitude Index), − miary statystyczne: odchylenie standardowe, wariancja, skoĹ›ność, kurtoza, współczynniki skali, ksztaĹ‚tu, statystyki wyĹźszego rzÄ™du,
31
? $ l M ^ =#
) = # # ] = # =# ) K* − wskaĹşniki minimalno-wariacyjne, − miary bazujÄ…ce na modelu, − alternatywne rozwiÄ…zania, np. wykorzystujÄ…ce sieci neuronowe, analizÄ™ w dziedzinie falek, czÄ™stotliwoĹ›ci, funkcji ortogonalnych, rozwiÄ…zania jÄ…drowe (kernelowe), fraktalne, analizy spektralne czy teĹź entropiÄ™, − czysto biznesowe wskaĹşniki, najczęściej wyraĹźone w jednostkach monetarnych zwane KPIs. W odpowiedzi na potrzebÄ™ doboru parametrĂłw regulatorĂłw PID opracowano szereg metod uĹ‚atwiajÄ…cych pracÄ™ automatykĂłw. Popularnymi metodami doboru parametrĂłw sÄ… [17, 25]: − strojenie rÄ™czne; − metody tabelaryczne (np. Ziegler-Nichols, Tyreus-Luyben); − strojenie automatyczne (np. relay tuning); − analiza czÄ™stotliwoĹ›ciowa; − optymalizacja wskaĹşnikĂłw jakoĹ›ci. Mimo Ĺźe wszystkie metody speĹ‚niajÄ… swoje zadanie, nie sÄ… pozbawione istotnych wad. Strojenie rÄ™czne zajmuje duĹźo czasu i nie daje Ĺźadnej gwarancji, Ĺźe uzyskane parametry sÄ… bliskie optymalnym. Metody automatyczne i ujÄ™te w tabelach reguĹ‚y wymagajÄ… specjalnego pobudzenia obiektu w trakcie procedury. Analiza czÄ™stotliwoĹ›ciowa z kolei pozwala zapewnić ukĹ‚adowi odpowiednie zapasy stabilnoĹ›ci, ale nie zapewnia bezpoĹ›redniego przeĹ‚oĹźenia na inne wskaĹşniki jakoĹ›ci, szczegĂłlnie te w dziedzinie czasu (czas sterowania, przeregulowania) ale i na wskaĹşniki caĹ‚kowe. Interpretacja zapasĂłw moduĹ‚u i fazy nie jest intuicyjna i trudna w przemysĹ‚owych zastosowaniach. Poza tym wymaga identyfikacji charakterystyki czÄ™stotliwoĹ›ciowej obiektu, co w przypadku wielu obiektĂłw (chemia procesowa, obiekty niestabilne) jest trudne lub wrÄ™cz niemoĹźliwe do uzyskania. Optymalizacja wskaĹşnikĂłw jakoĹ›ci, bÄ™dÄ…ca tematem niniejszej pracy, rĂłwnieĹź nie jest pozbawiona wad. Strojenie regulatora PID tÄ… metodÄ… wymaga przede wszystkim dokĹ‚adnego modelu procesu, ktĂłrego uzyskanie nie zawsze jest proste, a takĹźe wymaga niemaĹ‚ego nakĹ‚adu mocy obliczeniowej w samym procesie optymalizacji. Ponadto otrzymane w wyniku optymalizacji regulatory czÄ™sto dziaĹ‚ajÄ… zbyt agresywnie, Ĺźeby byĹ‚y bezpieczne dla urzÄ…dzeĹ„ wykonawczych wystÄ™pujÄ…cych w ukĹ‚adzie. Pomimo tych wad optymalizacja wskaĹşnikĂłw jakoĹ›ci umoĹźliwia dobĂłr parametrĂłw regulatora speĹ‚niajÄ…cych optimum jakoĹ›ci wedĹ‚ug okreĹ›lonego przez uĹźytkownika kryterium, a ponadto nie wymaga oddziaĹ‚ywania na obiekt w trakcie jego pracy. Co wiÄ™cej, wraz z rosnÄ…cÄ… mocÄ… obliczeniowÄ… sterownikĂłw programowalnych i komputerĂłw przemysĹ‚owych, a takĹźe coraz czÄ™stszym w przemyĹ›le archiwizowaniem przebiegĂłw zmiennych procesowych w bazach danych, wady tej metody powinny stawać siÄ™ coraz mniej dotkliwe w przyszĹ‚oĹ›ci [12]. RozwĂłj technik komputerowych i algorytmĂłw obliczeniowych jest motorem do odkrywania wĹ‚aĹ›ciwoĹ›ci algorytmĂłw poprzez ich symulacjÄ™ komputerowÄ…. Silne i powszechnie dostÄ™pne maszyny
pozwalają na wykonanie bardziej złoşonych obliczeń w krótkim czasie bez wymuszania stosowania technik pozwalających na analityczne uproszczenia. Idąc w tym nurcie autorzy niniejszej pracy postanowili zbadać wpływ wskaźników jakości opisanych na błędzie na parametry sygnału sterowania generowanego przez algorytm PID. Badania przeprowadzono w oparciu o symulacje układu regulacji i obiektu (zamknięty układ sterowania) – w kaşdym badanym przykładzie dokonano optymalizacji parametrów regulatora PID wykorzystując algorytm optymalizacji globalnej. Wybór ten ma na celu unikniecie potencjalnych rozwiązań lokalnych i często jest stosowany do takich zadań. W szczególności przebadano wskaźniki róşne od najczęściej stosowanego podejścia opartego na błędzie średniokwadratowym, którego popularność wynika właśnie z właściwości. Wskaźnik średniokwadratowy w swojej naturze jest zawsze dodatni i róşniczkowalny, dzięki czemu zastosowanie jego pozwala na uzyskanie pełnego analitycznego rozwiązania optymalnego lub teş ograniczeniu obliczeń w szeregu problemów [17]. Celem niniejszej pracy jest porównanie wybranych wskaźników jakości wykorzystywanych do strojenia regulatorów PID pod kątem właściwości sygnałów wyjściowych (sterowań) wyliczonych z algorytmu oraz próba znalezienia wskaźnika, który pozwala na uzyskanie bezpiecznych przebiegów sygnału sterującego dla wszystkich badanych klas obiektów. Treść została uporządkowana tak, aby moşliwie czytelnie zobrazować wykorzystane metody i algorytm oraz zaprezentować metodykę badań jakości sterowania. W oparciu o przeprowadzone badania sformułowana wnioski oraz wskazano moşliwe kierunki dalszych prac.
Z2 - # W badaniach zastosowano algorytm PID w strukturze rĂłwnolegĹ‚ej (rys. 1). W praktyce stosowana jest rzeczywista rĂłwnolegĹ‚a struktura regulatora PID, w ktĂłrej czĹ‚on róşniczkujÄ…cy jest rozszerzony o dziaĹ‚anie inercyjne ze staĹ‚Ä… czasowa Td/Îą. Regulator opisany jest wzorem (1): ⎛ ⎞ ⎜ 1 Td s âŽ&#x; âŽ&#x;E s + U s = K p ⎜1 + ⎜ Ti s Td s + 1 âŽ&#x; ⎜ âŽ&#x; Îą âŽ? âŽ
()
()
(1)
gdzie: U(s) – transformata sygnału sterowania; E(s) – transformata sygnału uchybu regulacji; Kp – współczynnik wzmocnienia; Ti – stała czasowa całkowania; Td – stała czasowa róşniczkowania; a – wzmocnienie dynamiczne regulatora, wartość jego zazwyczaj wybierana jest z zakresu 5–20, a wartość typowa wynosi 10.
Rys. 1. Struktura regulatora PID z rzeczywistym róşniczkowaniem Fig. 1. The structure of the PID controller with real differentiation
32
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
& ' ( $ ) *& * ) '
Rys. 2. Struktura układu regulacji z regulatorem PID Fig. 2. Structure of the control system with PID controller
Symulacje przeprowadzono w ten sposób, şe regulator PID podłączono do obiektu (opis badanych obiektów przedstawiono w rozdziale 2) w pętli ujemnego sprzęşenia zwrotnego (rys. 2). Celowo, co zostało wyjaśnione w rozdziale 4, w badaniach nie analizowano wpływu zakłóceń działających na obiekt. Dla takiego układu wyznaczono transmitancję zastępczą w sposób zgodny z rysunkiem 3.
Rys. 3. Schemat układu zamkniętego Fig. 3. Closed loop schema
Transmitancję układu zamkniętego moşna opisać wzorem (2).
()
L s =
() H (s ) + 1 H s
(2)
gdzie: Gr(s) – transmitancja regulator PID; Go(s) – transmitancja obiektu; L(s) – transmitancja zastępcza pętli regulacji; H(s) – transmitancja zastępcza regulatora PID i obiektu wyraşona jako: H(s) = Gr(s)¡Go(s). Uzyskany w ten sposób układ jest równowaşny początkowej postaci pętli regulacji i w takiej postaci został wykorzystywany w trakcie optymalizacji nastaw regulatorów. Uşycie zastępczej transmitancji pozwoliło na skrócenie czasu obliczeń w symulacjach, co w przypadku stosowania algorytmów optymalizacji, gdzie kaşda iteracja algorytmu optymalizacji wymaga wyliczenia wskaźnika jakości (i tym samym wykonania symulacji układu) miało istotne znaczenie. Celem zadania optymalizacji było wyznaczenie parametrów regulator PID, dla których wskaźnik jakości opisany na błędzie między wartością zadaną a uzyskanym przebiegiem osiągnie wartość minimalną. W symulacjach wykorzystano kilka typów wskaźników, opis wykorzystanych wskaźników przedstawiono w rozdziale 3. Dla kaşdego typu obiektu i kaşdego wskaźnika wykonano optymalizację nastaw algorytmu PID uşywając hybrydowego autorskiego rozwiązania polegającego na wstępnym siłowym przeszukaniu dziedziny rozwiązań po siatce, a następnie uruchomieniu funkcji środowiska MATLAB fmincon dla najlepszego, względem optymalizowanego wskaźnika, węzła z siatki. Kolejne kroki algorytmu w szczegółach przedstawiono w punktach ponişej: 1. Ręczne oszacowanie dziedziny parametrów regulatora PID, dla danego obiektu. 2. Określenie granulacji parametrów PID – gęstości siatki, dla której przeprowadzone zostaną symulacje. 3. Symulacja układu zamkniętego dla wszystkich kombinacji wartości parametrów regulatora PID (węzłów siatki) na
trajektorii w postaci skoku jednostkowego wartości zadanej o długości 500. 4. Wyznaczenie wartości uchybów dla kaşdego węzła i wyliczenie wskaźników jakości. 5. Wybór węzłów, dla których uzyskano najlepsze wartości względem badanych wskaźników jakości i uruchomienie, dla kaşdego wskaźnika algorytmu lokalnego fmincon z rozwiązaniem startowym równym wartościom parametrów (nastaw regulatora PID) z węzła. Opisany algorytm zastosowano dla kaşdego badanego obiektu. W ten sposób uzyskano optymalne zestawy parametrów strojeniowych algorytmu PID niezaleşnie dla kaşdego badanego obiektu i kaşdego badanego wskaźnika jakości. Znalezione rozwiązania optymalne przebadano następnie pod kątem właściwości sygnału sterującego, jak i uzyskanych przebiegów sygnałów procesowych.
v2 7 ) * W celu uwiarygodnienia i sprawdzenia uniwersalności sformułowanych wniosków, badania przeprowadzono na obiektach o róşnych właściwościach. Obiekty zostały wybrane ze względu na swoja uniwersalność i moşliwość modelowania przy ich pomocy wielu rzeczywistych procesów przemysłowych, co na przykład widać w klasycznych przykładach dla układów regulacji typu PID [2]. Podczas badań wykorzystano pięć obiektów (z czterech klas – obiekt oscylacyjny badany był z opóźnieniem i bez opóźnienia) i przebadano sześć róşnych wskaźników jakości, według których dokonano optymalizacji parametrów algorytmu PID.
v2Y2 7 W artykule opisano rezultaty z badaĹ„ obiektĂłw opisanych rĂłwnaniami (3–7): − obiekt inercyjny I rzÄ™du
()
G s =
1 10s + 1
(3)
− obiekt inercyjny II rzÄ™du z opóźnieniem
()
G s =
1 e −10s 100s 2 + 40s + 1
(4)
− obiekt oscylacyjny i obiekt oscylacyjny z opóźnieniem
()
1 4s 2 + 0,8s + 1
(5)
1 e −10s 4s 2 + 1,2s + 1
(6)
G s =
()
G s = −
obiekt nieminimalnofazowy
()
G s =
−10s + 1 400s 2 + 40s + 1
(7)
33
? $ l M ^ =#
) = # # ] = # =# ) K*
v2Z2 ) *
− wartość minimalna Umin, − maksymalny przyrost ΔUmax, − wariancja sterowania ĎƒU, − Ĺ›redni przyrost sterowania A(ΔU), − mediana przyrostĂłw sterowania M(ΔU).
Wykorzystano sześć wskaźników jakości opisanych równaniami (8–13). Wybór łączy w sobie wskaźniki klasyczne z odpornymi. Wynika to z faktu, şe charakter przebiegów często cechuje się właściwościami niegaussowskimi [6, 7], a takich przypadkach wskaźniki odporne jak IAE, LMS czy teş ADM powinny być lepiej dopasowane. ISE (Integral of Squared Error) – całka z kwadratu uchybu regulacji na całym horyzoncie czasowym eksperymentu [17], wskaźnik ten opisany jest wzorem (8):
âˆŤ
ISE =
tk
t0
e 2dt
(8)
IAE (Integral of Absolut Error) – całka z modułu uchybu regulacji na całym horyzoncie czasowym eksperymentu [18], wskaźnik ten opisany jest wzorem (9):
IAE =
âˆŤ
tk
t0
e dt
(9)
LMS (Least Median of Squares) – mediana kwadratów uchybu regulacji [15], wskaźnik ten opisany jest wzorem (10):
( )
LMS = median e 2
Rys. 4. Trajektoria zadana do badani sygnału sterującego Fig. 4. Trajectory given to test control signal
(10)
ADM (Average Distance from Median) – średnia odległość próbki od mediany [24], wskaźnik ten opisany jest wzorem (11):
(
( ))
ADM = x1 − M x
W tym miejscu naleşy zauwaşyć, şe wszelkie analizy prowadzone są dla działania nadąşnego pętli regulacji w odpowiedzi na zmiany wartości zadanej. Niemniej nic nie stoi na przeszkodzie w rozszerzeniu interpretacji wyników na zadania tłumienia zakłóceń. Wskaźniki oceny jakości sterowania wyznaczane są dla uchybu regulacji. Z punktu widzenia dynamicznej jakości regulacji, źródło owych zmian nie jest najwaşniejsze. Regulator musi jak najszybciej je zredukować. Tym samym metodologia oceny układów sterowania wydaje się być analogicznie stosowalna. Takowe badania obejmujące zmiany zarówno wartości zadanej jak i zakłóceń moşna równieş znaleźć w literaturze [5]. Warto równieş zauwaşyć, şe w prowadzonych badaniach nie występują şadne nieliniowości czy teş ograniczenia. Pętla regulacji jest całkowicie liniowa. Celem pracy jest pokazanie, şe nawet przy takich jakşe idealistycznych załoşeniach istnieje sensowność i celowość stosowania wskaźników odpornych i nie-Gaussowskich. Wprowadzenie do pętli nieliniowości, jakimi są ograniczenia sygnały sterującego moşe tylko bardziej przyczynić się do sensowności stosowania wskaźników nietradycyjnych jakimi są wskaźniki odporne.
(11)
ITSE (Integral of Time-weighted Squared Error) – modyfikacja wskaźnika ISE polegająca na całkowaniu iloczynu kwadratu uchybu i czasu, dla którego ten uchyb wystąpił [17], wskaźnik opisany jest wzorem (12):
ITSE =
âˆŤ (t − t )e dt tk
2
(12)
0
t0
ITAE (Integral of Time-weighted Absolut Error) – analogiczna do ITSE modyfikacja wskaźnika IAE polegająca na całkowaniu iloczynu wartości bezwzględnej uchybu i czasu, dla którego ten uchyb wystąpił [17], wskaźnik ten opisany jest wzorem (13):
ITAE =
âˆŤ (t − t ) e dt tk
(13)
0
t0
x2 - * Â
' # (
X2
W wyniku procedury opisanej w podrozdziale „Metody i algorytmyâ€? uzyskane zostaĹ‚y nastawy optymalnych pod kÄ…tem wybranych wskaĹşnikĂłw jakoĹ›ci regulatorĂłw PID. Jako, Ĺźe do optymalizacji wymagana byĹ‚a znajomość jedynie wartoĹ›ci zadanej i wyjĹ›cia obiektu, dla uproszczenia obliczeĹ„ wykorzystywana byĹ‚a transmitancja pÄ™tli zamkniÄ™tej, a takĹźe stosunkowo krĂłtka trajektoria zadana o dĹ‚ugoĹ›ci 500 krokĂłw (czas dyskretyzacji wynosiĹ‚ 1 s). Jako trajektorii zadanej uĹźyto sygnaĹ‚ skoku jednostkowego, zmiana wartoĹ›ci zadanej nastÄ™powaĹ‚a w kroku jedenastym. Po znalezieniu szukanych parametrĂłw regulatora PID, jego sygnaĹ‚ sterujÄ…cy jest badany po przeprowadzeniu symulacji na duĹźo dĹ‚uĹźszej i bardziej zróşnicowanej trajektorii zadanej (rys. 4). Uzyskane przebiegi zostaĹ‚y zapisane, a nastÄ™pnie obliczono nastÄ™pujÄ…ce miary zmiennoĹ›ci sygnaĹ‚u sterujÄ…cego: − wartość maksymalna Umax,
34
P
O
M
I
A
R
Y
•
A
U
T
O
M
Dla kaĹźdego obiektu i kaĹźdego wskaĹşnika przeprowadzono symulacje, wyniki przedstawiono poniĹźej w postaci: − optymalnych nastaw regulatora PID; − wybranych fragmentĂłw przebiegĂłw odpowiedzi ukĹ‚adu na skok jednostkowy i sygnaĹ‚u sterowania; − wartoĹ›ci wskaĹşnika dla trajektorii zróşnicowanej. Ze wzglÄ™du na czytelność zdecydowano o przedstawieniu tylko wybranych fragmentĂłw przebiegĂłw odpowiedzi ukĹ‚adu na skok jednostkowy i sygnaĹ‚u sterowania.
X2Y2 7 ' # 0 W tabeli 1 przedstawiono parametry regulatorów PID dla poszczególnych wskaźników jakości. Na rysunkach 5 i 6 zaprezentowano przebiegi sygnałów sterującego i sterowanego dla przebadanych wskaźników jakości. W tabeli 2 przedstawiono wartości miar zmienności sygnału sterującego. A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
& ' ( $ ) *& * ) '
Rys. 5. Przebiegi sygnałów w pętli regulacji dla regulatora optymalizowanego pod kątem wskaźników: ISE, IAE, ADM, ITSE, ITAE Fig. 5. Process and control signals for the controller optimized for indicators: ISE, IAE, ADM, ITSE, ITAE
Rys. 6. Przebiegi sygnałów w pętli regulacji dla regulatora optymalizowanego pod kątem wskaźnika LMS Fig. 6. Process and control signals for the controller optimized for LMS indicator
Tabela 1. Parametry regulatorĂłw PID uzyskane dla obiektu inercyjnego Table 1. Parameters of PID controllers obtained for the inertial object
wartości zadanej zajmuje mu około 20 s. Bardzo dobra jakość regulacji większości badanych regulatorów PID wynika z charakterystyki badanego obiektu. Pojedyncza inercja bez opóźnienia jest stosunkowo prostym obiektem do sterowania. Jeśli chodzi o parametry sygnałów sterujących, regulator strojony wskaźnikiem LMS ponownie znacząco odbiega od pozostałych. Obliczone przez ten regulator maksymalne i minimalne wartości sterowania są ponad 5-krotnie mniejsze od wszystkich pozostałych wskaźników. Średni przyrost oraz wariancja sterowania równieş są dla tego regulatora kilkukrotnie nişsze niş w przypadku pozostałych. Jednak, co ciekawe, mediana przyrostów sterowania dla regulatora strojonego LMS nie jest juş tak niska i regulatory strojone wskaźnikami ISE, a tym bardziej ADM osiągają pod tym kątem lepsze wyniki. Moşe to wynikać z wkładu członu róşniczkującego, jako şe regulatory z niską stałą czasową róşniczkowania Td mają najnişsze mediany przyrostów sterowania. Najlepsze parametry sygnału sterującego poza regulatorem strojonym wskaźnikiem LMS uzyskał PID wyznaczony w oparciu o ITSE, a następnie kolejno IAE oraz ADM. Co ciekawe, mimo şe wskaźnik IAE pozwolił uzyskać lepsze rezultaty niş ISE, to w odpowiadających im wskaźnikach waşonych po czasie lepsze przebiegi sterowania od ITAE pozwolił uzyskać wskaźnik ITSE.
Kp
Ti
Td
ISE
98,8132
9,9583
0,0321
IAE
85,3976
10,0281
0,0280
LMS
5,0000
1,0000
0,3000
ADM
98,9823
10,0123
0,0123
ITSE
93,0511
9,9860
0,0095
ITAE
76,0477
10,0824
0,0816
Przebiegi sygnału wyjściowego procesu niewiele się miedzy sobą róşnią, za wyjątkiem regulatora optymalizowanego wskaźnikiem LMS. Podczas gdy w pozostałych układach regulacji sygnał wyjściowy obiektu szybko osiąga wartość zadaną i nie powoduje przeregulowań, regulator optymalizowany LMS ma przeregulowanie sięgające 28% a ustabilizowanie procesu na
Tabela 2. Właściwości sygnału sterującego uzyskanego dla obiektu inercyjnego Table 2. Properties of the control signal obtained for the inertial object
Umin
Umax
ΔUmax
A(ΔU)
M(ΔU)
ĎƒU
ISE
–103,6210
130,5322
130,5322
0,2420
0,0068
9,3827
IAE
–86,6448
109,3089
109,3089
0,2073
0,0164
6,6052
LMS
–15,1914
20,0000
20,0000
0,0448
0,0081
0,3465
ADM
–88,1116
111,1571
111,1571
0,2047
0,0038
6,8266
ITSE
–109,6694
138,1026
138,1026
0,2596
0,0085
10,4963
ITAE
–80,7120
101,8910
101,8910
0,1825
7,695×10–15
5,7485
35
? $ l M ^ =#
) = # # ] = # =# ) K*
X2Z2 7 ++ 0 ' )
Na rysunku 7 zaprezentowano przebiegi sygnałów sterującego i sterowanego dla przebadanych wskaźników jakości. W tabeli 4 przedstawiono wartości miar zmienności sygnału sterującego. W porównaniu z obiektem inercyjnym badany obiekt o charakterystyce inercji drugiego rzędu z opóźnieniem jest duşo trudniejszym procesem do wysterowania. Znajduje to odzwierciedlenie w uzyskanych sygnałach wyjściowych z obiektu sterowania, wszystkie regulatory poza PID optymalizowanym przy pomocy wskaźnika ITAE powodują wyraźne przeregulowania i oscylacje sygnału wyjściowego. Z kolei regulator uzyskany w wyniku optymalizacji wskaźnika ITAE powoduje aperiodyczną odpowiedz układu (nie występuje przeregulowanie), a przy tym sprawia, şe proces osiąga wartość zadaną w podobnym czasie, co w przypadku pozostałych regulatorów. Generalnie zapewnia zatem zdecydowanie najlepszą jakość regulacji dla tego obiektu. Z punktu widzenia uzyskanych sygnałów sterujących regulator optymalizowany wskaźnikiem ITAE ponownie okazał się najlepszy. Generowany przez niego sygnał sterujący ma najmniejsze wartości maksymalne i minimalne, najmniejszy maksymalny i średni przyrost sterowania, a takşe najnişszą wariancję. Pod względem mediany przyrostów sterowania najlepszy okazał się regulator uzyskany w wyniku optymalizacji wskaźnika LMS, a lepsze od ITAE okazały się równieş regulatory optymalizowane wskaźnikami ADM oraz IAE. Stosunkowo wysoka mediana przyrostów sterowania w przypadku PID uzyskanego z optymalizacji wskaźnika ITAE moşe wynikać z bardzo niskiej w porównaniu z innymi regulatorami stałej czasowej róşniczkowania Td.
W tabeli 3 przedstawiono parametry regulatorów PID dla poszczególnych wskaźników jakości. Tabela 3. Parametry regulatorów PID uzyskane dla obiektu dwuinercyjnego z opóźnieniem Table 3. Parameters of PID controllers obtained for the two-inertial object with delay
Kp
Ti
Td
ISE
2,8434
36,1037
7,1016
IAE
2,6669
42,9624
5,2216
LMS
3,0000
27,0000
4,0000
ADM
2,6675
43,0032
5,2188
ITSE
2,8314
38,8992
6,3860
ITAE
1,6368
38,3431
0,0330
X2v2 7 ' ) W tabeli 5 przedstawiono parametry regulatorów PID dla poszczególnych wskaźników jakości. Na rysunkach 8 i 9 zaprezentowano przebiegi sygnałów sterującego i sterowanego dla przebadanych wskaźników jakości. W tabeli 6 przedstawiono wartości miar zmienności sygnału sterującego. Podobnie jak w przypadku obiektu inercyjnego pod względem przebiegu sygnału wyjściowego najbardziej od pozostałych odróşnia się regulator otrzymany w wyniku optymalizacji wskaźnika LMS. Tak samo jak pozostałe regulatory nie powoduje on istotnego przeregulowania, jednak ustabilizowanie procesu w poblişu wartości zadanej zajmuje mu około 15 s, podczas gdy pozostałe regulatory robią to prawie natychmiastowo. Ponadto w trakcie narastania sygnału wyjściowego
Rys. 7. Przebiegi sygnałów w pętli regulacji dla regulatora optymalizowanego pod kątem badanych wskaźników Fig. 7. Process and control signals for the controller optimized for the investigated indicators
Tabela 4. Właściwości sygnału sterującego uzyskanego dla obiektu dwuinercyjnego z opóźnieniem Table 4. Properties of the control signal obtained for the two-inertial object with delay
36
Umin
Umax
ΔUmax
A(ΔU)
M(ΔU)
ĎƒU
ISE
–24,2219
31,2774
31,2774
0,0656
3,5211×10–7
0,7611
IAE
–22,6687
29,3359
29,3359
0,0560
4,3536×10–8
0,6247
LMS
–25,6000
33,0000
33,0000
0,0624
7,8840×10–11
0,7540
ADM
–22,6740
29,3425
29,3425
0,0560
3,9704×10–8
0,6248
ITSE
–24,1163
31,1454
31,1454
0,0630
8,3855×10–8
0,7276
ITAE
–13,6038
18,0048
18,0048
0,0323
6,2593×10–8
0,3028
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
& ' ( $ ) *& * ) '
procesu moşna zaobserwować wyraźne oscylacje związane z charakterystyką obiektu. Nastawy regulatorów uzyskanych w wyniku optymalizacji wskaźników jakości innych niş LMS są bardzo zblişone, zatem podobieństwo otrzymanych przebiegów sygnałów wyjściowych sterowanego procesu oraz sygnałów sterujących nie jest zaskakujące. Mimo bardziej złoşonej charakterystyki obiektu regulacji, regulatory doprowadzają proces do wartości zadanej błyskawicznie nie powodując przy tym istotnych przeregulowań. Okupione jest to jednak ogromnymi skokami sygnału sterującego. W porównaniu z regulatorem optymalizowanym wskaźnikiem LMS pozostałe regulatory wygenerowały sygnały sterujące o wartościach minimalnych i maksymalnych wyşszych ponad 75-krotnie. Podobnie wygląda kwestia maksymalnego oraz średniego przyrostu sterowania. Mediana przyrostów sterowania i wariancja sterowania w przypadku regulatora optymalizowanego LMS okazały się mniejsze od uzyskanych przez pozostałe regulatory nawet o 3 rzędy wielkości.
Tabela 5. Parametry regulatorów PID uzyskane dla obiektu oscylacyjnego bez opóźnienia Table 5. Parameters of PID controllers obtained for the oscillation object without delay
Kp
Ti
Td
ISE
78,1200
1,9433
1,3898
IAE
77,8801
1,7351
1,4447
LMS
1
2
2,4000
ADM
77,8740
1,7619
1,4650
ITSE
78,0604
1,8772
1,4088
ITAE
77,8361
1,6690
1,4278
Rys. 8. Przebiegi sygnałów w pętli regulacji dla regulatora optymalizowanego pod kątem wskaźników ISE, IAE, ADM, ITSE, ITAE Fig. 8. Process and control signals for the controller optimized for indicators: ISE, IAE, ADM, ITSE, ITAE
Rys. 9. Przebiegi sygnałów w pętli regulacji dla regulatora optymalizowanego pod kątem wskaźnika LMS Fig. 9. Process and control signals for the controller optimized for LMS indicator
Tabela 6. Własności sygnału sterującego uzyskanego dla obiektu oscylacyjnego bez opóźnienia Table 6. Properties of the control signal obtained for the oscillation object without delay
Umin
Umax
ΔUmax
A(ΔU)
M(ΔU)
ĎƒU
ISE
–687,1566
859,3200
865,5620
2,1691
0,7140
404,6738
IAE
–684,4369
856,6811
864,0537
1,8326
0,2468
401,8422
LMS
–7,9997
11,0000
11,1462
0,0213
1,3572×10–4
0,1761
ADM
–684,4185
856,6140
864,7042
1,7240
0,1747
401,7528
ITSE
–686,3749
858,6644
865,3759
2,2886
0,8187
403,7944
ITAE
–684,2061
856,1971
862,8163
1,8832
0,2805
401,4695
37
? $ l M ^ =#
) = # # ] = # =# ) K*
X2x2 7 ' )
Na rysunku 10 zaprezentowano przebiegi sygnałów sterującego i sterowanego dla przebadanych wskaźników jakości. W tabeli 8 przedstawiono wartości miar zmienności sygnału sterującego. Wprowadzenie opóźnienia do obiektu oscylacyjnego sprawiło, şe uzyskane nastawy regulatorów róşnią się o wiele wyraźniej niş w przypadku obiektu bez opóźnień. Obiekt oscylacyjny z duşym opóźnieniem jest wyraźnie trudniejszy do ustabilizowania. Jedynym regulatorem, przy którym obiekt nie wpada w drgania i nie powoduje przeregulowania jest PID uzyskany w wyniku optymalizacji wskaźnika LMS. Ten regulator zapewnia asymptotyczny przebieg wyjścia procesu, a takşe stosunkowo szybko (choć wolniej niş pozostałe regulatory) doprowadza do osiągnięcia wartości zadanej. Poza uzyskaniem dobrej jakości regulacji procesu, generowany przez regulator optymalizowany wskaźnikiem LMS sygnał sterujący jest zdecydowanie najwolniej zmienny ze wszystkich uzyskanych dla tego obiektu. Najmniejsza i największa wartość sterowania są ponad dwukrotnie mniejsze niş w przypadku najlepszego z pozostałych regulatorów, podobnie jak maksymalny i średni przyrost sterowania. Mediana przyrostów sterowania jest 8 rzędów wielkości mniejsza niş pozostałe, a wariancja równieş jest najnişsza, choć róşnica w stosunku do pozostałych regulatorów jest mniej wyraźna. Poza regulatorem LMS, warto wyróşnić regulator ITAE, który wywołuje najmniejsze spośród zaobserwowanych przeregulowań, a takşe szybko osiąga wartość zadaną. Ponadto pod względem zmienności sterowania ustępuje tylko regulatorowi optymalizowanemu wskaźnikiem LMS.
Ponişej w tabeli przedstawiono parametry regulatorów PID dla poszczególnych wskaźników jakości. Tabela 7. Parametry regulatorów PID uzyskane dla obiektu oscylacyjnego z opóźnieniem Table 7. Parameters of PID controllers obtained for the oscillation object with delay
Kp
Ti
Td
ISE
0,3916
4,4275
3,1632
IAE
0,3783
4,6034
2,4319
LMS
0,1500
2,9000
0,8000
ADM
0,3783
4,6033
2,4319
ITSE
0,3887
4,5180
2,8746
ITAE
0,3533
4,5604
2,0362
X2X2 7
% Ponişej w tabeli 9 przedstawiono parametry regulatorów PID dla poszczególnych wskaźników jakości. Na rysunku 11 zaprezentowano przebiegi sygnałów sterującego i sterowanego dla przebadanych wskaźników jakości. W tabeli 10 przedstawiono wartości miar zmienności sygnału sterującego. Dla obiektu nieminimalnofazowego uzyskano najbardziej zróşnicowane nastawy regulatorów ze wszystkich obserwowanych obiektów. Najmniej agresywny sygnał sterujący generował ponownie regulator uzyskany na drodze optymalizacji wskaźnika LMS. Ten regulator powodował najmniejszy maksymalny przyrost sterowania, najmniejszy średni przyrost sterowania, najmniejsze skrajne wartości sygnału sterującego oraz najmniejszą wariancję. Natomiast mediana przyrostów sterowania dla tego regulatora okazała się największa ze wszystkich dla tego obiektu. Niestety łagodny przebieg sygnału sterującego w tym
Rys. 10. Przebiegi sygnałów w pętli regulacji dla regulatora optymalizowanego pod kątem badanych wskaźników Fig. 10. Process and control signals for the controller optimized for the investigated indicators
Tabela 8. Własności sygnału sterującego uzyskanego dla obiektu oscylacyjnego z opóźnieniem Table 8. Properties of the control signal obtained for the oscillation object with delay
38
Umin
Umax
ΔUmax
A(ΔU)
M(ΔU)
ĎƒU
ISE
–2,6461
4,3076
4,3076
0,0118
3,9023 10–7
0,1234
IAE
–2,5296
4,1613
4,1613
0,0106
0,0013
0,1221
LMS
–0,5195
1,6500
1,6500
0,0035
1,1102 10–15
0,1144
ADM
–2,5296
4,1613
4,1613
0,0106
0,0013
0,1221
ITSE
–2,6203
4,2757
4,2757
0,0109
0,0003
0,1229
ITAE
–2,3081
3,8863
3,8863
0,0084
0,0002
0,1210
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
& ' ( $ ) *& * ) '
Tabela 9. Parametry regulatorĂłw PID uzyskane dla obiektu nieminimalnofazowego Table 9. Parameters of PID controllers obtained for the non-minimum phase object
Kp
Ti
Td
ISE
1,8986
39,9297
9,5907
IAE
2,4198
39,1037
9,5016
LMS
1,0791
18,8611
6,1131
ADM
3
39
8,6000
ITSE
2,1786
39,3102
9,8247
ITAE
2,4407
39,1113
9,3505
przypadku wiązał się z wystąpieniem przeregulowania, oscylacji i wolnego osiągnięcia wartości zadanej przez regulator. Drugi pod względem zmienności sterowania okazał się regulator ISE, trzeci ITSE, czwarty IAE, piąty ITAE, a najbardziej agresywny sygnał sterujący generował regulator optymalizowany wskaźnikiem ADM. Parametry sygnałów sterowania są ściśle powiązane ze wzmocnieniami statycznymi Kp uzyskanych regulatorów, a ich odzwierciedlenie widać bardzo wyraźnie w uzyskanych odpowiedziach obiektu sterowania. Moşna by uznać, şe im większe jest wzmocnienie regulatora PID i bardziej agresywne generowane przez niego sterowanie, tym krótszy czas regulacji, a większe przeregulowanie, jednak regulator optymalizowany wskaźnikiem LMS jest wyjątkiem od tej prawidłowości. Warto zwrócić uwagę na fakt, şe dotychczas w kaşdym badanym obiekcie sterowania regulator optymalizowany w oparciu o wskaźnik ISE uzyskiwał bardziej agresywne sterowanie niş ten, który był optymalizowany przy zastosowaniu wskaźnika IAE, natomiast w przypadku obiektu minimalnofazowego jest dokładnie odwrotnie. Podobnie minimalizacja wskaźnika ITAE dawała do tej pory łagodniejsze przebiegi sygnału sterującego w porównaniu z ITSE. Inną godną zastanowienia obserwacją jest fakt, şe optymalizacja przy wykorzystaniu wskaźnika ADM po raz pierwszy dała wyraźnie odbiegający od innych zestaw parametrów PID, podczas gdy do tej pory były one bardzo zblişone do regulatorów optymalizowanych wskaźnikiem IAE.
V2 & U *
Rys. 11. Przebiegi sygnałów w pętli regulacji dla regulatora optymalizowanego pod kątem badanych wskaźników Fig. 11. Process and control signals for the controller optimized for the investigated indicators
Celem pracy było porównanie przebiegów czasowych w układach regulacji PID o róşnych nastawach regulatorów oraz wartości wybranych kryteriów jakości regulacji obliczanych na podstawie tych przebiegów i cel ten zrealizowano. Niestety uzyskane wyniki są często właściwe poszczególnym klasom obiektów i dokonanie na ich podstawie syntezy uniwersalnych stwierdzeń prawdziwych dla dowolnego układu regulacji jest trudne. Mimo şe zaleşności między parametrami sygnałów wyjściowych regulatorów zmieniają się dla róşnych właściwości sterowanych procesów, niektóre trendy da się zaobserwować dla wszystkich rozpatrywanych obiektów regulacji. Symulacje na obiekcie inercyjnym pokazały, şe dla takiego prostego procesu, wszystkie regulatory optymalizowane całkowymi wskaźnikami sprowadzały obiekt na wartość zadaną w sposób natychmiastowy, jednak odbywało się to kosztem bardzo wysokich i gwałtownych przyrostów sygnału sterują-
Tabela 10. Właściwości sygnału sterującego uzyskanego dla obiektu nieminimalnofazowego Table 10. Properties of the control signal obtained for the non-minimum phase object
Umin
Umax
ΔUmax
A(ΔU)
M(ΔU)
ĎƒU
ISE
–15,9077
20,8846
20,8846
0,0380
2,7852 10–10
0,4912
IAE
–20,4942
26,6178
26,6178
0,0522
3,0342 10–10
0,8531
LMS
–8,6963
11,8701
11,8701
0,0218
2,1098 10–4
0,2204
ADM
–25,6000
33
33,0000
0,0695
1,9457 10–9
1,4119
ITSE
–18,3717
23,9646
23,9646
0,0457
7,1368 10–12
0,6782
ITAE
–20,6782
26,8477
26,8477
0,0523
0,0002
0,8552
39
? $ l M ^ =#
) = # # ] = # =# ) K*
cego w momencie zmiany wartości zadanej. Podobnie zachował się PID powstały w wyniku minimalizacji wskaźnika ADM. Zblişone rezultaty uzyskano dla obiektu oscylacyjnego. Jedyna róşnicą były marginalne oscylacje pojawiające się na wyjściu sterowanego obiektu, oraz jeszcze bardziej agresywne przebiegi sygnału sterującego. Z kolei regulatory o nastawach dobieranych w oparciu o wskaźnik LMS zarówno w przypadku obiektu inercyjnego, jak i oscylacyjnego zapewniały sygnał sterujący o wielokrotnie nişszych niş w przypadku pozostałych regulatorów wartościach krańcowych, zdecydowanie wolniejszych przyrostach oraz nişszej wariancji. Niestety w obu przypadkach wiązało się to ze spadkiem jakości regulacji – wolniejszym osiąganiem wartości zadanej, a w przypadku obiektu oscylacyjnego takşe wyraźnym przeregulowaniem. Dopiero obiekty o wolniejszej odpowiedzi na zmiany sygnału wejściowego, czyli takie w których występują opóźnienia lub nieminimalnfazowość sprawiły, şe róşnice między badanymi wskaźnikami jakości stały się wyraźniejsze. W przypadku rozpatrywanego obiektu nieminimalnofazowego oraz obiektu oscylacyjnego z opóźnieniem najbardziej łagodny przebieg sygnału sterującego ponownie zapewniły regulatory PID o nastawach wyznaczonych na drodze minimalizacji wskaźnika LMS. Natomiast dla obiektu o charakterystyce inercji drugiego rzędu z opóźnieniem najwolniej zmienny sygnał sterujący wygenerował regulator optymalizowany pod kątem wskaźnika ITAE, podczas gdy sterowanie uzyskane przez regulator strojony względem LMS było najbardziej agresywne ze wszystkich zaobserwowanych. Prawdopodobnie wynika to z niskiej stałej czasowej całkowania Ti, która w związku z duşym opóźnieniem w sterowanym procesie powoduje wysoki i gwałtowny przyrost sterowania po zmianie wartości zadanej. Moşliwe, şe uzyskane wyniki mają związek z monotonicznością odpowiedzi skokowej odpowiednich obiektów i wskaźnik LMS daje najłagodniejsze przebiegi sterowania regulatorów tylko, gdy obiekt regulacji ma niewielkie opóźnienie lub jego odpowiedz skokowa jest niemonotoniczna. Związek miedzy regulatorami optymalizowanym pod kątem ISE oraz IAE jest stały dla wszystkich badanych obiektów, które są minimalnofazowe. PID uzyskane w wyniku optymalizacji ISE są dla tych obiektów zawsze bardziej agresywne w porównaniu z IAE. Najmniejsze róşnice są widoczne dla obiektu oscylacyjnego, a najbardziej dla obiektu dwuinercyjnego z opóźnieniem. Analogicznie regulatory uzyskane w wyniku optymalizacji ITSE generują bardziej agresywny sygnał sterujący niş PID strojone względem ITAE. Dla wszystkich minimalnofazowych obiektów ITAE daje łagodniejsze sterowanie niş IAE, natomiast miedzy ITSE a ISE nie ma şadnej stałej zaleşności. Zaleşności między wskaźnikami ISE a IAE oraz między ITSE a ITAE są dokładnie odwrotne w przypadku badanego obiektu nieminimalnofazowego, Zachowanie to moşe wynikać z właściwości uchybu występującego zaraz po zmianie wartości zadanej – w przypadku obiektu nieminimalnofazowego uchyb wzrasta powyşej wartości 1, gdy dla wszystkich innych obiektów uchyb jest zawsze co najwyşej równy wartości 1. Niniejsza praca nie wyczerpuje tematu, jest szereg zagadnień, które wpływają na charakter sygnałów sterowania i byłyby warte badania. Jednym z zagadnień jest uwzględnienie odporności optymalizowanych regulatorów i zapasów stabilności. Moşliwe byłoby zmodyfikowanie wskaźników jakości w sposób taki, şeby w procesie optymalizacji uwzględnić powyşsze aspekty i zbadać jaki jest wpływ ich na uzyskane przebiegi sygnałów sterowania i jakość regulacji.
40
P
O
M
I
A
R
Y
•
A
U
T
O
M
Ciekawą propozycją byłoby uşycie wskaźników wielokryterialnych. Tym sposobem moşliwe byłoby wykorzystanie na przykład wskaźnika LMS lub podobnych do ograniczania nagłych zmian sygnału sterowania, a bardziej agresywnego wskaźnika, takiego jak ISE lub ISE w celu zwiększenia szybkości i dokładności regulacji procesu. Modyfikując wagi poszczególnych wskaźników projektant układów sterowania korzystający z takiego rozwiązania miałby moşliwość dostosowania nastaw regulatorów miedzy bardziej agresywnymi, a generującymi łagodny przebieg sterowania w zaleşności od potrzeb.
~ # % 1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
A
T
Ă…strĂśm K.J., Computer control of a paper machine – an application of linear stochastic control theory, “IBM Journal of Research and Developmentâ€?, Vol. 11, No. 4, 1967, 389–405, DOI: 10.1147/rd.114.0389. Ă…strĂśm K.J., Hägglund T., Benchmark Systems for PID Control, “IFAC Proceedings Volumesâ€?, Vol. 33, No. 4, 2000, 165–166, DOI: 10.1016/S1474-6670(17)38238-1. Bauer M., Horch A., Xie L., Jelali M., Thornhill N., The current state of control loop performance monitoring – a survey of application in industry, “Journal of Process Controlâ€?, Vol. 38, 2016, 1–10, DOI: 10.1016/j.jprocont.2015.11.002. Choudhury M.A.A.S., Shah S.L., Thornhill N.F., Diagnosis of Process Nonlinearities and Valve Stiction. Advances in Industrial Control. Springer Berlin Heidelberg, 2008. Doma ski P.D., Non-Gaussian properties of the real industrial control error in SISO loops, [in:] Proceedings of the 19th International Conference on System Theory, Control and Computing, Cheile Gradistei, Romania, 2015, 877– 882, DOI: 10.1109/ICSTCC.2015.7321405. Doma ski P.D., Non-Gaussian statistical measures of control performance, “Control and Cyberneticsâ€?, Vol. 46, No. 3, 2017, 259–290. Doma ski P.D., Golonka S., Jankowski R., Kalbarczyk P., Moszowski B., Control rehabilitation impact on production nefficiency of ammonia synthesis installation , “Industrial & Engineering Chemistry Researchâ€?, Vol. 55, 2016, 10366– 10376, DOI: 10.1021/acs.iecr.6b02907. Gao X., Yang F., Shang C., Huang D., A review of control loop monitoring and diagnosis: Prospects of controller maintenance in big data era, “Chinese Journal of Chemical Engineeringâ€?, Vol. 24, No. 8, 2016, 952–962, DOI: 10.1016/j.cjche.2016.05.039. Harris T., Assessment of closed loop performance, “The Canadian Journal of Chemical Engineeringâ€?, Vol. 67, 1989, 856–861, DOI: 10.1002/cjce.5450670519. Horch A., Isaksson A.J., A modified index for control performance assessment, [in:] Proceedings of the 1998 American Control Conference, 1998, 3430–3434, DOI: 10.1109/ACC.1998.703231. Jelali M., An overview of control performance assessment technology and industrial applications. “Control Engineering Practiceâ€?, Vol. 14, No. 5, 2006, 441–466, DOI: 10.1016/j.conengprac.2005.11.005. Jelali M., Control Performance Management in Industrial Automation: Assessment, Diagnosis and Improvement of Control Loop Performance, Springer-Verlag, London 2013. Ordys A., Uduehi D., Johnson M.A., Process Control Performance Assessment – From Theory to Implementation. Springer-Verlag, London 2007.
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
& ' ( $ ) *& * ) '
14. Paulonis M.A., Cox J.W., A practical approach for largescale controller performance assessment, diagnosis, and improvement, “Journal of Process Controlâ€?, Vol. 13, No. 2, 2003, 155–168, DOI: 10.1016/S0959-1524(02)00018-5. 15. Rousseeuw P.J., Leroy A.M., Robust Regression and Outlier Detection, John Wiley & Sons, Inc., New York, NY, USA, 1987. 16. Schäfer J., Cinar A., Multivariable MPC system performance assessment, monitoring, and diagnosis. “Journal of Process Controlâ€?, Vol. 14, No. 2, 2004, 113–129, DOI: 10.1016/j.jprocont.2003.07.003. 17. Seborg D.E., Mellichamp D.A., Edgar T.F., Doyle F.J., Process dynamics and control, Wiley, 2010. 18. Shinskey F.G., How good are our controllers in absolute performance and robustness? “Measurement and Controlâ€?, Vol. 23, No. 4, 1990, 114–121, DOI: 10.1177/002029409002300402. 19. Smuts J.F., Hussey A., Requirements for successfully implementing and sustaining advanced control applications, [in:] Proceedings of the 54th ISA POWID Symposium, Charlotte, North Carolina, USA, 89–105. 20. Srinivasan B., Spinner T., Rengaswamy R., Control loop performance assessment using detrended fluctuation analysis (DFA), “Automaticaâ€?, Vol. 48, No. 7, 2012, 1359–1363, 2012, DOI: 10.1016/j.automatica.2012.04.003. 21. Starr K.D., Petersen H., Bauer M., Control loop performance monitoring – ABB’s experience over two decades, “IFAC-PapersOnLineâ€?, Vol. 49, No. 7, 2016, 526–532, 11th IFAC Symposium on Dynamics and Control of Process Systems Including Biosystems DYCOPS-CAB 2016, Trondheim, Norway, DOI: 10.1016/j.ifacol.2016.07.396. 22. Tatjewski P., Sterowanie zaawansowane procesĂłw przemysĹ‚owych. Struktury i algorytmy, Akademicka Oficyna Wydawnicza EXIT, Warszawa 2016. 23. Tolfo F., A methodology to assess the economic returns of advanced control projects. 1983 American Control Conference, IEEE, 1141–1146, DOI: 10.23919/ACC.1983.4788287.
24. Verboven S., Hubert M., LIBRA: a MATLAB library for robust analysis, “Chemometrics and Intelligent Laboratory Systems�, Vol. 75, No. 2, 2005, 127–136, DOI: 10.1016/j.chemolab.2004.06.003. 25. Veronesi M., Visioli A., An industrial application of a performance assessment and retuning technique for PI controllers, “ISA Transactions�, Vol. 49, No. 2, 2010, 244–248, DOI: 10.1016/j.isatra.2009.11.008. 26. DeVries W., Wu S., Evaluation of process control effectiveness and diagnosis of variation in paper basis weight via multivariate time-series analysis. “IEEE Transactions on Automatic Control�, Vol. 23, No. 4, 1978, 702–708, DOI: 10.1109/TAC.1978.1101828. 27. Zhuo H., Research of performance assessment and monitoring for multivariate model predictive control system. [in:] 4th International Conference on Computer Science & Education, 2009, 509–514, DOI: 10.1109/ICCSE.2009.5228377.
4 ] # % 4 % 0 # ( 046 The goal of this thesis is to assess and compare PID controllers with parameters determined by minimizing select performance indices paying special attention to the attributes of their output signals. Analysis is based on simulations performed using MATLAB for four controlled processes classes and six types of indicators. The methods and algorithms used have been presented in detail. Thesis presents results of the simulations, optimized PID parameters, plots of process signals in examined control systems, calculated control signal attributes and formulated based on experiments observations and conclusions. Keywords` 046 & & # & # % & (
41
? $ l M ^ =#
) = # # ] = # =# ) K*
" 4 $
" ) $ >
F 0 ' %
' %
8% 4 8 4 # = > % ? = . 4 # 0 > 8 F S 6 868 ; < F % T %< = (< ( %
8F 8 G F > % ? = . 4 # 0 = > S % 44 = F S % & = (
! " " $ # 5" 5 >
% ' % 4 @ & % 4 8 4 # > % ? . 4 # = 0 > 0 S = % S _ = ( F S S = %< = % = %
42
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
Pomiary Automatyka Robotyka, ISSN 1427-9126, R. 23, Nr 1/2019, 43–46, DOI: 10.14313/PAR_231/43
B % % % % F F @ F % 5( $
V B % & G @ *& !)=-+E >
4 $
I A % Y "8 4?>4 N 0 4 8 0 < 0480& 8 R )!)& !)=-KJ >
) ? $ 0 A & > % 4 @ S% & $ . )& *J=!!* "
4 W artykule przedstawiono nową koncepcję modelu doskonalenia jakości w przedsiębiorstwie produkcyjnym branşy biomedycznej. Zaproponowane rozwiązanie uwzględnia zarówno specyficzne wymagania branşy biomedycznej, jak równieş integrację procesu doskonalenia jakości z rozwiązaniami zarządzania zautomatyzowanymi systemami produkcji. Zaproponowany system doskonalenia jakości obejmuje ponadto oryginalną metodę monitorowania odpadów z produkcji oraz elementów wadliwych (czerwone oko), która nie tylko zapewnia poprawę efektywności produkcji, ale równieş ułatwia działania pro-ekologiczne i umoşliwia zmniejszenie obciąşenia środowiska. Zaproponowany model został zweryfikowany w zakładzie produkującym soczewki wewnątrzgałkowe, prowadząc do wzrostu produktywności o ponad 20%. ( ` % & F % & I
1. Wprowadzenie Od początku XXI wieku w krajach wysoko rozwiniętych moşna zaobserwować intensyfikację procesów społecznych związanych z problemami starzejącego się społeczeństwa. Równocześnie, starzenie się społeczeństwa prowadzi do gwałtownego wzrostu wydatków na słuşbę zdrowia, w szczególności na okulistykę. Przykładowo wydatki na operacje zaćmy w Polsce w latach 2015–2017 wzrosły trzykrotnie [1]. W rezultacie starzenia się społeczeństwa, gwałtownie wzrosło zapotrzebowanie na zaawansowane produkty biomedyczne, takie jak na przykład soczewki wewnątrzgałkowe. Produkcja wyrobów biomedycznych jest jednym z najbardziej technologicznie zaawansowanych obszarów inşynierii produkcji. W odróşnieniu na przykład od produkcji w branşy motoryzacyjnej, produkcja biomedyczna jest przedmiotem ścisłych regulacji prawnych, określających zarówno nadzór nad półproduktami, sam proces produkcji, transport, magazynowanie oraz utylizację zuşytych wyrobów biomedycznych.
. ' `
G & ' . ( % !* ! )!*+ & % % )E ! )!*+
Tak ścisłe regulacje prawne powodują konieczność rozwoju specjalizowanych metod doskonalenia jakości produkcji wyrobów biomedycznych. W doskonaleniu procesów produkcji biomedycznej naleşy uwzględnić wymagania prawne i normatywne, które czynią tę branşę jedną z najbardziej konserwatywnych branş produkcyjnych we współczesnej gospodarce. Równocześnie poprawa jakości procesu produkcji przedsiębiorstwa branşy biomedycznej prowadzi zarówno do znaczącego wzrostu jego rentowności, jak równieş umoşliwia lepsze zaspokojenie potrzeb odbiorców produktu – zarówno lekarzy jak i pacjentów. Mając na uwadze istotę poruszanego zagadnienia w artykule przedstawiono nową koncepcję modelu doskonalenia jakości, który uwzględni potrzeby branşy biomedycznej, a w szczególności spełni rygorystyczne wymogi prawne.
Z2 / ' '
* Model doskonalenia jakości musi być ściśle związany z systemem jakości przedsiębiorstwa, obejmującym nie tylko proces produkcyjny, lecz takşe inne kluczowe procesy w nim realizowane. Na rysunku 1 przedstawiono uogólniony schemat procesów w zakładzie produkcyjnym. W odniesieniu do tak zdefiniowanego samego procesu produkcji, jak i innych procesów zachodzących w zakładzie produkcyjnym zaproponowano model doskonalenia systemu jakości. Model ten przedstawiono schematycznie na rysunku 2.
!
43
#
#
$ $ ! $ )
Wymagania prawne Badania satysfakcji klienta zapobiegawcze
Analiza danych procesowych
Proces produkcji Zasoby -k adrowe -i nfrastrukturalne
Produkt
intelektualnej
Wymagania
Informacje o produkcie
Proces technologiczny Rys. 1. Schemat procesów zachodzących w firmie powiązanych z modelem doskonalenia jakości produktu Fig. 1. Schematic block diagram of processes in production plant, connected to the quality increase model
Zaproponowany model doskonalenia produkcji w branşy biomedycznej bazuje na mapie procesu [2] wykonanej zgodnie ze standardami przemysłowymi. Następnie proces w sposób ciągły analizowany jest z wykorzystaniem metod subiektywnych i obiektywnych. W zakresie metod subiektywnych, w praktyce stosuje się sformalizowany system sugestii [3] oraz metodę „5 S� [4], której nazwa pochodzi od japońskich słów: Seiri (sortowanie), Seiton (systematyzacja), Seiso (sprzątanie), Seiketsu (standaryzacja), Shitsuke (samodyscyplina). Liczba poszczególnych działań, oznaczanych literą „S� w praktyce japońskich firm, skąd pochodzi metoda, mogła być róşna, ze względu na uwarunkowania panujące w organizacji [5]. Na przykład, w procesie produkcyjnym firmy Toyota nie uwzględniano analizy w zakresie samodyscypliny [6]. W odniesieniu do metod obiektywnych w systemach produkcji biomedycznej stosuje się analizę OEE (ang. Overall Equipment Effectiveness) oraz monitorowanie kluczowych parametrów wydajności procesu [7] KPI (ang. Key Performance Indicators). W analizie OEE przyjmuje się, şe [8]:
Podstawowe załoşenie metody czerwonego oka to organizacja spotkania produkcyjnego, podczas którego omawiane są wszystkie wadliwe komponenty, półwyroby, które nie spełniły norm jakościowych. Podczas tego spotkania omawiane są wady, które ujawniły się podczas procesu wytwarzania. Na bazie analizy odpadów oraz wyników analizy obiektywnej i subiektywnej, właściciele obszarów – kierownicy i liderzy podejmują działania korygujące i naprawcze. Naleşy podkreślić, şe ze względu na specyfikę branşy biomedycznej działania korygujące i naprawcze muszą mieć charakter interdyscyplinarnych projektów, w które zaangaşowani są nie tylko technolodzy, lecz takşe specjaliści zorientowani na zagadnienia certyfikacji procesu produkcyjnego wyrobów oraz spełnienia wymagań normatywnych. Coraz częściej w działaniach korygujących i naprawczych wiodącą rolę odgrywają rzecznicy patentowi prowadzący w sposób ciągły badania stanu techniki w odniesieniu do procesów produkcyjnych.
3. Wyniki wdroşenia zaproponowanego modelu Zaproponowany model doskonalenia produkcji w branşy biomedycznej został wdroşony w 2013 r. w warszawskim zakładzie produkcyjnym firmy Valeant Med. Wśród korzyści wprowadzenia zaproponowanych rozwiązań metody naleşy podkreślić poprawę kultury technicznej i kultury organizacji pracy, poniewaş metoda ta angaşuje wszystkich pracowników do procesu rozwiązywania problemu. Ponadto pracownicy produkcji nie próbują rozwiązywać problemu sami, lecz uzyskują znaczące wsparcie działów wspomagających. Z ekonomicznego punktu widzenia, dzięki wprowadzeniu proponowanego modelu doskonalenia produkcji doprowadziło do wzrostu produktywności rozumianej jako udział czasu pracy przeznaczonego na wytwarzanie w całkowitym czasie eksploatacji linii produkcyjnej. Wzrost produktywności w zakładzie przedstawiono ilościowo w tabeli 1.
OEE = (B/A) Ă— (D/C) Ă— (F/E) Ă— 100% przy czym kolejne czynniki tego iloczynu to: − dostÄ™pność (B/A), gdzie B to rzeczywisty czas pracy danej jednostki technologicznej, zaĹ› A to oczekiwany czas jej pracy; − efektywność dziaĹ‚ania (D/C), gdzie C to docelowa wydajność, zaĹ› D jest wydajnoĹ›ciÄ… aktualnÄ…; − jakość okreĹ›lona przez (F/E), gdzie E to caĹ‚kowita liczba wyprodukowanych sztuk, zaĹ› F liczba sztuk zgodnych z wymaganiami. Informacje pozyskane zarĂłwno z analizy subiektywnej jak i obiektywnej stanowiÄ… dane wejĹ›ciowe dla metody czerwonego oka. Metoda ta bazuje na kole Deminga [9], jednak w szczegĂłlny sposĂłb koncentruje siÄ™ na analizie odpadĂłw i wadliwych produktĂłw oraz komponentĂłw.
44
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
* $ " ) ( @#
Mapa procesu
Analiza subiektywna
Analiza obiektywna
5S
OEE Golem Analiza KPI
Metoda Czerwonego Oka
zgodnie z wymaganiami prawnymi i normatywnymi Rys. 2. Schemat proponowanego modelu doskonalenia produkcji w branĹźy biomedycznej Fig. 2. Schematic block diagram of proposed model for quality increase in biomedical production plant
Tabela 1. Produktywność warszawskiego zakładu firmy Valeant Med w latach 2012–2017 Table 1. Productiveness of Warsaw division of Valeant Med Company in years 2012–2017
Rok
2012
2013
2014
2015
2016
2017
Produktywność
59,7%
74,4%
79,3%
82,9%
82,8%
80,8%
x2 &
Jak wykazano w artykule, w rezultacie starzenia się społeczeństwa, gwałtownie wzrosło zapotrzebowanie na zaawansowane produkty biomedyczne. W odpowiedzi na to zapotrzebowanie korzystne będzie podjęcie działań zmierzających do wzrostu produktywności w branşy biomedycznej. Zaproponowany model doskonalenia produkcji w branşy biomedycznej uwzględnia jej specyficzne wymagania prawne i normatywne. Bazując zarówno na analizie subiektywnej, jak i obiektywnej, oraz na analizie odpadów z procesu produkcji, zaproponowany model umoşliwia prawidłową organizację działań korygujących i naprawczych. Działania te muszą mieć charakter interdyscyplinarnych projektów, w które zaangaşowani są technolodzy, specjaliści od zagadnień certyfikacji procesu produkcyjnego i wyrobów oraz rzecznicy patentowi. Poprawność działania proponowanego modelu doskonalenia produkcji w branşy biomedycznej została zweryfikowana praktycznie w warszawskim oddziale firmy Valeant Med. W wyniku wdroşenia zaproponowanego modelu, w latach 2012–2017 wskaźnik produktywności w oddziale wzrósł o ponad 20%.
~ # % 1. Grzybowski A., Maciejewski A., Koszty społeczne odroczonego terminu zaćmy z jednoczesnym wszczepianiem soczewki, OphthaTherapy, Vol. 3, Nr 1(9), 2016, 53–58.
2. Grajewski P., Procesowe zarzÄ…dzanie organizacjÄ…, Polskie Wydawnictwo Ekonomiczne, Warszawa 2012. 3. Grycuk A., Dekier Ĺ ., Programy sugestii pracowniczych. DoĹ›wiadczenia polskich przedsiÄ™biorstw, Stowarzyszenia Lean Management Polska, WrocĹ‚aw 2014. 4. Jaca C., Viles E., Paipa-Galeano L., Santos J., Mateo R., Learning 5S principles from Japanese best practitioners: case studies of five manufacturing companies, “International Journal of Production Researchâ€?, Vol. 52 No. 15, 2014, 4574–4586, DOI: 10.1080/00207543.2013.878481. 5. Gapp R., Fisher R., Kobayashi K., Implementing 5S within a Japanese context: an integrated management system, “Management Decisionâ€? Vol. 46, No. 4, 2008, 565–579, DOI: 10.1108/00251740810865067. 6. Ohno T., Taichi Ohno’s workplace management. Special 100th birthday edition, McGraw Hill, 2013. 7. Torkko M., Katajavuori N., Linna A., Juppo A.M., The Utilization of Quality KPIs in the Pharmaceutical Industry, “Journal of Pharmaceutical Innovationâ€?, Vol. 8, No. 3, 2014, 175–182, DOI: 10.1007/s12247-014-9184-3. 8. Sonmez V., Testik M.C., Testik O.M., Overall equipment effectiveness when production speeds and stoppage durations are uncertain, “The International Journal of Advanced Manufacturing Technologyâ€?, Vol. 95, No. 1–4, 2018, 121–130, DOI: 10.1007/s00170-017-1170-8. 9. Deming W.E., Out of the Crisis, MIT Press, 1986.
45
#
#
$ $ ! $ )
. B % # 6 # 0 % 0 A % 0 % 0 Paper presents the new concept of quality improvement model for biomedical production plant. Proposed solutions take into consideration both specific requirements of biomedical production as well as integration of quality assessment process with automated production control systems. Moreover, proposed quality improvement model covers original method of production waste monitoring (red eye method), which provide not only increase of production efficiency but also opens new possibilities of ecology-oriented actions. Proposed model was tested in intraocular lenses production process leading to the increase of productiveness exceeding 20%. Keywords` % & F % Z & Z
" 5( $
! " " 4 $
% F '
'
8F > % B 0 = > )!!! $% )! 6 % % % F = & F = % $% % % % V B % % S ( !! !!! S ( P 0 X ( % R % F = 0 <% N < ( = " S 4 F ; % S
> A % = Y "8 4?>4 N 0 = ( 4 8 0 < 0480 0 % = S ) ; G F % # = ( > ( 0 % . # 0 = & B G L)!*JN)!*QM 0 % S ( ; G F = *- % < L F = < M ( B ; % 0 ( $ ( 4 ( G < L)!*EN)!*KM " % = % < < # ( < %= F & < % # ( = % <% ( F @ & F 0 " % % < < 0 ( A % T LX A GM ( & < % F
) ? $ ( % ' F % 0 > % 4 @
S% 0 A = 8 FS%^ < = < & ( # ( < F %S # % F % T % # S T ( & S% ( S& ( # % & F % T # ( ( ( & S% S " = ( < F % ( S ( = % ( F F % T # ( = & < F % < F S # % & < ( & # F 0 ( . S% 0 = % S& 4 # B # (& % B ( ? ( (& 4??? ` ( 0 # 4??? 0 % . ( B (
46
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
Pomiary Automatyka Robotyka, ISSN 1427-9126, R. 23, Nr 1/2019, 47–58, DOI: 10.14313/PAR_231/47
% % I * 0 % % Zygmunt Lech Warsza I A % Y "8 4?>4 N 0 4 8 0 < 0480& 8 R )!)& !)=-KJ >
Jacek Puchalski ; < S% B & ? )& !!=!!* >
Streszczenie: W dwuczęściowej pracy przedstawia się zmodyfikowaną wersję wektorowej metody oceny niepewności wieloparametrowych pomiarów pośrednich. Jest to rozszerzenie metody podanej w Suplemencie 2 do Przewodnika wyznaczania niepewności pomiarów. Nowością zaproponowaną w niniejszej pracy jest uwzględnianie wpływu skorelowania rozkładów mierzonych wartości elementów multimenzurandu na niepewności składowe typu A lub/oraz typu B mierzonych pośrednio parametrów menzurandu wyjściowego. Omówienia dokonano na przykładzie pomiarów menzurandu dwuparametrowego 2D o skorelowanych parametrach mierzonych. W części 1. wyznaczono wzory ogólne dla macierzy kowariancji, niepewności i współczynnika korelacji menzurandu wyjściowego oraz dla kilku charakterystycznych przypadków szczególnych. Na wykresach zilustrowano zaleşności współczynnika korelacji wielkości wyjściowych w funkcji udziałów składowych typu B w niepewnościach wielkości wejściowych. W części 2. będą omówione przykłady estymacji niepewności i współczynnika korelacji dla kilku rodzajów funkcji przetwarzania wyników pomiarów oraz sformułuje się wnioski ogólne. ( ` & % )6& ( & % < 8 A& % <
1. Wprowadzenie Wyniki pomiarów podawane są zwykle w postaci dwu estymat, które opisują wartości oraz oceny dokładności kaşdej z wielkości badanego obiektu, czyli parametry menzurandu. Oba estymatory wyznacza się z danych próbki pomiarowej, czyli zbioru wartości obserwacji pomiarowych, które powtarzano wielokrotnie w tych samych lub podobnych warunkach wpływających. Podczas tworzenia oceny dokładności pomiarów wieloparametrowych uwzględnia się związki deterministyczne i statystyczne (korelację) między mierzonymi wielkościami. Po korekcie znanych wpływów o charakterze deterministycznym otrzymuje się rozrzut wartości obserwacji pomiarowych traktowany juş jako
. ' `
( P > & *+ J'( . ( % *+ !) )!*+ & % % )K ! )!*+
!
losowy i opisywany rozkładem prawdopodobieństwa. W celu ujednolicenia oceny dokładności wyników pomiarów, siedem organizacji międzynarodowych wspólnie opracowało i pierwszy raz opublikowało w 1993 r. specjalny dokument – Przewodnik wyraşania niepewności pomiarów (ang. Evaluation of measurement data — Guide to the expression of uncertainty in measurement), znany pod angielskim akronimem GUM [1]. Jego treść udoskonalano kilkakrotnie. Ostatnia wersja ukazała się w 2008 r. [1]. Zakres stosowania przewodnika GUM jest rozszerzany przez kolejne Suplementy [2–5]. Ciągle ukazują się prace o usprawnieniu korzystania i udoskonaleniu tych przepisów, w tym i przyczynki z udziałem autorów [8, 10–20]. Poprzednio do oceny dokładności pomiarów stosowano błędy pomiarowe z podziałem na systematyczne i przypadkowe. Wskutek trudności z wyznaczeniem w praktyce wartości rzeczywistej, w przewodniku GUM wprowadzono nowe pojęcie nazwane niepewnością pomiarów (ang. measurement uncertainty). Jest to szerokość przedziału, a dla pomiarów wieloparametrowych – opis granic tzw. obszaru pokrycia, w którym z określonym prawdopodobieństwem moşe wystąpić estymata wartości menzurandu, którą otrzymano po korekcji i przetworzeniu surowych wyników pomiarów. Podstawą oceny niepewności pomiarów jest wyznaczenie jej składowych typu A i typu B, oznaczonych jako uA i uB. Niepewność całkowita u jest zdefiniowana w GUM
47
W ) ^
) ) & ; ^H U& &&& jako odchylenie standardowe rozkĹ‚adu wypadkowego w postaci splotu dwĂłch niezaleĹźnych statystycznie rozkĹ‚adĂłw skĹ‚adowych. Pierwszy z nich opisuje statystycznie rozrzut otrzymanych eksperymentalnie wartoĹ›ci obserwacji pomiarowych. Drugi jest rozkĹ‚adem hipotetycznym randomizujÄ…cym przypuszczalne zmiany wynikĂłw róşnego pochodzenia, w tym o charakterze systematycznym, nieznane co do wartoĹ›ci w trakcie wykonywanych pomiarĂłw. Zmiany te mogÄ… pojawiać siÄ™ losowo w dĹ‚ugim okresie uĹźytkowania mierzonego obiektu, przyrzÄ…du lub systemu pomiarowego i w róşnych dopuszczalnych dla nich warunkach otoczenia, innych niĹź przy kalibracji przyrzÄ…du i przeprowadzaniu pomiarĂłw kontrolnych. Natomiast zwykle nie zmieniajÄ… siÄ™ w trakcie pomiarĂłw. Funkcje gÄ™stoĹ›ci prawdopodobieĹ„stwa rozkĹ‚adu wypadkowego i rozkĹ‚adĂłw skĹ‚adowych speĹ‚niajÄ… warunek p(x) = p(xA) ∗ p(xB). Niepewność caĹ‚kowita, czyli odchylenie standardowe u rozkĹ‚adu wypadkowego jest sumÄ… geometrycznÄ… niepewnoĹ›ci skĹ‚adowych uA i uB, tj.:
u = uA2 + uB2
maksymalne, czy teş przez stosowane nadal jeszcze błędy graniczne. Na podstawie ich wartości, przy załoşeniu równomiernego rozkładu, wyznacza się odchylenie standardowe jako składową niepewności typu B przyrządu. Właściwości metrologiczne aparatury pomiarowej decydują o wartościach zmierzonych nią parametrów wielu urządzeń i procesów spełniających odpowiedzialne funkcje. Stosowane w przemyśle i innych działach gospodarki systemy jakości, obok starannej kalibracji aparatury pomiarowej wymagają teş okresowej kontroli metrologicznej, jak i kontroli podstawowych parametrów badanych procesów i urządzeń w trakcie okresu ich eksploatacji [6, 7, 20]. W przeprowadzanych w praktyce eksperymentach pomiarowych bada się teş parametry skojarzonych ze sobą menzurandów jedno- i wieloparametrowych występujących w tym samym, lub róşnych obiektach mierzonych. Tylko niekiedy wykonuje się je jednym, a zwykle wieloma przyrządami, w tych samych lub róşnych warunkach otoczenia. Otrzymywany rozrzut wartości powtarzanych obserwacji w pomiarach wieloparametrowych jest wywołany zarówno wspólnymi jak i róşnymi przyczynami. Wielowariantowość sytuacji powoduje, şe całkowite i tzw. rozszerzone, o określonym prawdopodobieństwie niepewności wielkości mierzonych jako elementów multimenzurandu, zaleşą nie tylko od ich składowych typu A oraz typu B, ale są teş ze sobą skojarzone wskutek wzajemnego skorelowania między tymi wielkościami. Podobnie jak niepewność typu A, skorelowanie moşe przejawiać się między wartościami obserwacji uzyskiwanymi w trakcie wykonywania danego eksperymentu pomiarowego, jak teş tak jak niepewność typu B powiązać ze sobą wyniki pomiarów róşnych eksperymentów wykonywanych w innych warunkach. Zagadnień tych nie obejmują w pełni rekomendacje szacowania niepewności pomiarów zawarte w GUM i w jego Suplemencie 2 [3] oraz w innych przepisach międzynarodowych, np. [6, 7]. Celem tej dwuczęściowej publikacji jest prezentacja rozszerzonej metody wyznaczania niepewności w pośrednich pomiarach wieloparametrowych, gdy wielkości mierzone na wejściu są skorelowane. Uwzględnia się w niej wyodrębnienie wpływu skorelowania wielkości mierzonych na poszczególne składowe ich niepewności typu A oraz/lub B. Takiej metody dotąd nie opisywano w literaturze. W pracy przeprowadzi się analizę wpływu róşnych przypadków skorelowania wielkości mierzonych jako elementów multimenzurandu na macierze kowariancji dla składowych niepewności typu A oraz B oraz na wypadkową macierz kowariancji. Zilustruje się to kilkoma przykładami pomiarów dwuparametrowych opisanych funkcją liniową i prostymi funkcjami nieliniowymi. Wyznaczy się teş odpowiednie wzory. Jest to istotne szczególnie, gdy wyniki wyjściowe pomiarów wieloparametrowych będą potem uşywane wspólnie.
(1)
Następnie dla znanej funkcji opisującej gęstość prawdopodobieństwa rozkładu, np. funkcji Gaussa dla rozkładu normalnego, na podstawie odchylenia standardowego i współczynnika rozszerzenia kP [1], wyznacza się analitycznie niepewność rozszerzoną U o określonym prawdopodobieństwie pokrycia P, np. 95% lub 99%. Dla innych otrzymanych eksperymentalnie rozkładów danych pomiarowych, niepewność rozszerzoną U wyznacza się teş metodą numeryczną Monte Carlo wg Suplementu 1 [2] do przewodnika GUM. Niepewność składową uA, zaleşną od rozkładu obserwacji pomiarowych, wyznacza się znanymi metodami statystycznymi, jak dla rozkładu normalnego [1] lub innych rozkładów [8, 21]. Natomiast niepewność uB szacuje się heurystycznie na podstawie subiektywnej wiedzy o przewidywanych zakresach i rozkładach wielkości wpływających i funkcjach ich oddziaływań. Niepewność uB randomizuje moşliwe oddziaływania wielu takich wielkości, których bieşące wartości nie są znane i brakuje danych, aby obliczyć dla nich poprawki. Wpływy te są nieusuwalne z wartości obserwacji pomiarowych. Zacharow przy wyznaczaniu niepewności pośrednich pomiarów jednoparametrowych rozpatrywał szacowanie heurystyczne skorelowania wielkości wpływających na niepewność typu B, które nazywał „logicznym� [9]. Dotyczy ono sytuacji, gdy skorelowania nie moşna wyznaczyć bezpośrednio z danych pomiarowych. Dorozhovetz i Warsza uzasadnili [10, 11], şe w warunkach pracy o znanym ograniczonym zakresie zmian wielkości wpływających, moşna przyjąć istotnie mniejszą wartość składowej uB, niş podana po kalibracji przyrządu dla pełnego zakresu dopuszczalnych zmian warunków jego pracy. Międzynarodowy Komitet ds. Przewodników w Metrologii ICGM proponował ostatnio opracowanie nowej wersji przewodnika GUM 2, opartej na prawdopodobieństwie warunkowym wg podejścia Bayesa [8, Dodatek 2], ale wyznaczania niepewności typu B jeszcze nie omawiano. Niepewność pomiaru MU (ang. Measurement Uncertainty) ma teş ogromne zastosowanie jako podstawowy element oceny zgodności produktów z wymaganiami [5, 7, 8]. Umoşliwia porównywanie wyników róşnych testów, sprawdzanie przekraczania ograniczeń lub spełniania wymagań tolerancji w produkcji i przez produkty. Dzięki moşliwości porównania wyników kalibracji z wymaganiami stanowi teş podstawę metrologicznej akceptacji wyposaşenia pomiarowego. Laboratoria akredytowane zgodnie z ISO/IEC 17025 [6, 7] muszą obliczyć MU dla kaşdej metody badawczej stosowanej w przyznanym im zakresie akredytacji. Istnieje jednakşe pewien rozdźwięk między opisem dokładności pomiarów za pomocą niepewności i opisem dokładności przyrządów i urządzeń pomiarowych przez dopuszczalne błędy
48
P
O
M
I
A
R
Y
•
A
U
T
O
M
Z2
- ' ' #
' *
' ' * Przetwarzanie sygnałów pomiarowych w pomiarach wieloparametrowych i wyznaczanie wyjściowej macierzy kowariancji UY wyjaśniają schematy (rys. 1). Wielkości badane pośrednio, zwane teş obserwablami [13–15], tworzą wieloparametrowy menzurand wyjściowy Y. Korelacja przejawia się zarówno w trakcie danego eksperymentu, jak i między pomiarami wykonywanymi w róşnych warunkach wpływających. Wartości elementów multimenzurandu Y i jego macierz kowariancji UY otrzymuje się pośrednio z przetwarzania wyników bezpośrednio mierzonych wielkości wejściowych, czyli elementów xi multimenzurandu X. Wielkości mierzone w danym eksperymencie pomiarowym mają niepewności typu A, oznaczane przez uA, które wynikają z rozrzutu wartości obserwacji w próbce. Skorelowanie danych pomiarowych tych wielkości będzie opisywane współczynnikami korelacji z literą A w indeksie dolnym. Współczynniki korelacji opisujące dodatkowe skoreA
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
< =)# ? #
Wektorowa metoda wg Suplementu 2 do GUM
Rys. 1. Schematy przetwarzania sygnałów mierzonych i wyznaczania macierzy kowariancji w rozszerzonej metodzie propagacji niepewności dla pośrednich pomiarów wieloparametrowych: UAB – macierz kowariancji wielkości wejściowych X o niepewnościach typu A lub/i typu B uwzględniająca skorelowania, X i Y multimenzurandy: wejściowy i wyjściowy; UX i UY – ich macierze kowariancji; Z = G(Y) – następne przetwarzanie multimenzurandu wyjściowego Y Fig. 1. Schemes of the signal values processing and covariance matrix estimation in the extended method of uncertainty propagation for multiparameter measurements: UAB covariance matrix of input of quantities X with uncertainty components of type A or/and B and correlation coefficients, X and Y – input and output quantities; UX i UY their covariance matrixes, Z = G(Y) – next conversion of measurand Y
lowanie między wielkościami mierzonymi w dwu eksperymentach o róşnych warunkach wpływających na niepewność wypadkową typu B, będą miały w indeksie literę B. Opisuje się to wspólnie w macierzy UAB. Przy jej uşyciu wyznacza się macierz kowariancji UX dla multimenzurandu wejściowego X, a z niej, dla znanej funkcji przetwarzania Y = F(X) – macierz kowariancji UY. Macierz ta uwzględnia skorelowania wpływające na niepewności typu A lub/i typu B na wejściu. Drugim wariantem opisującym wpływ obu rodzajów skorelowania jest przedstawienie macierzy kowariancji jako sumy dwu macierzy, np. UX = UXA + UXB. Omówi się to w ostatnim rozdziale tej pracy. W pomiarach wieloparametrowych wyjściowy m-wymiarowy menzurand Y zaleşy pośrednio od mierzonego na wejściu n-wymiarowego menzurandu X wg ogólnej funkcji o postaci Y = F(X)
Do wyznaczania niepewności w pomiarach pośrednich wielowymiarowych menzurandów wykorzystuje się wektorowe prawo propagacji niepewności podane w Suplemencie 2 [3] do przewodnika GUM. Propagację niepewności opisuje się jako zaleşność między macierzami kowariancji wielkości X i Y występującą przy linearyzacji funkcjonału F ze wzoru (2) przez pochodne, tj.: (3)
gdzie macierz czuĹ‚oĹ›ci S oraz macierze kowariancji UX, UY wielkoĹ›ci wejĹ›ciowych X i wyjĹ›ciowych Y, opisane sÄ… nastÄ™pujÄ…cymi wzorami: ⎥ ∂y1 ⎢ ⎢ ∂x1 ⎢ S =⎢! ⎢ ⎢ ∂ym ⎢ ∂x ⎣ 1
UX
∂y1 ⎤ ⎼ ∂x n ⎼ ⎼ ! ! ⎼ ⎼ ∂ym ⎼ ! ∂x n ⎼⎌ !
⎥ ux21 ! Ď x 1nux 1uxn ⎤ ⎢ ⎼ ⎼ =⎢ ! ! ! ⎢ ⎼ 2 âŽ˘Ď u u ⎼ u ! xn ⎣ xn 1 xn x 1 ⎌
(3c)
Wartości elementów wektora Y mogą być teş dalej wspólnie przetwarzane wg innej funkcji: Z = G (Y )
(4)
Wówczas dla wyznaczanych kolejno wartości elementów, nowego multimenzurandu Z = [z1, ‌, zm]T i elementów jego macierzy kowariancji UZ z niepewnościami uz1, ‌, uzm stosuje się nowe równanie wektorowe o takiej samej postaci jak (3), opisujące propagację niepewności, tj.:
(2)
gdzie: Y i X – wektory o elementach, które są wartościami wielkości: menzurandu wyjściowego i wejściowego.
UY = SUXST
UY
⎥ uy21 ! Ď y 1muy 1uym ⎤ ⎢ ⎼ ⎼ =⎢ ! ! ! ⎢ ⎼ 2 âŽ˘Ď u u ⎼ u ! y 1 m ym y 1 ym ⎣ ⎌
(3a)
(3b)
UZ = SGUYSGT
(5)
gdzie SG jest nową macierzą czułości dla funkcji G, utworzoną analogicznie jak S. Wzory dla pomiarów pośrednich menzurandów 3D i kilka przykładów stosowania metody wektorowej w elektrycznych układach pomiarowych DC i AC autorzy przedstawili w pracach [16–20]. Wyniki pomiarów pośrednich wieloparametrowych wyznacza się z parametrów bezpośrednio mierzonego menzurandu X, otrzymanych po przetworzeniu wielokrotnie powtarzanych obserwacji traktowanych w ogólnym przypadku jako losowe i uzyskanych w określonych warunkach otoczenia. W tym celu dokonuje się obróbki statystycznej otrzymanych surowych danych pomiarowych oraz koryguje się te niepoşądane znane wpływy warunków otoczenia przez poprawki. Warunki te mogą mieć charakter stały albo w czasie pomiarów zmieniać się deterministycznie w znany sposób lub losowo. Wprowadzając poprawki eliminuje się znane błędy systematyczne. Pozostałe przyczyny niedokładności szacuje się przez niepewności typu B. Następnie wyznacza się niepewności typu A i typu B oraz niepewności wypadkowe u wszystkich parametrów menzurandu X ze wzoru (1) i tworzy się jego macierz kowariancji UX, w której uwzględnia się skorelowanie parametrów. Wartości parametrów badanych menzurandów Y oraz Z uzyskuje się pośrednio po przetworzeniu wg wzorów (2)–(5). Istnieje wiele moşliwości skorelowania wielkości będących elementami multimenzurandu. Moşna wyróşnić wpływy skorelowania dotyczące osobno kaşdego typu niepewności A i B tych
49
W ) ^
) ) & ; ^H U& &&&
W ogólnym przypadku wyniki pomiarów wykonywanych dla otrzymania wartości estymatorów x1, x2 będą miały róşne wartości par niepewności składowych u1A, u2A oraz u1B, u2B. Jeśli eksperymenty te często odbywają się w podobnych lub nawet w jednakowych warunkach, to wartości wielkości mierzonych mogą być skorelowane i powiązane współczynnikiem korelacji 1 ³ rB ≠0. Podobnie przy niejednakowych, ale częściowo wspólnych przyczynach losowego rozrzutu obserwacji pomiarowych dwu róşnych badań tego samego obiektu, w ogólnym przypadku mogą pojawić się teş niepewności uA i współczynnik skorelowania 1 ³ rA ≠0. Jeśli korelacja występuje między zbiorami odchyleń od wartości estymatorów wielkości, to standardowe odchylenia ich rozkładów są niepewnościami jednego typu – A lub B. Związki te ilustruje rysunek 2.
wielkości. Rozrzut wartości obserwacji pomiarowych zaleşy nie tylko od warunków otoczenia systemu pomiarowego i obiektu badanego, ale i od zmian parametrów ich elementów wewnętrznych. Próbki pomiarowe tego samego badanego obiektu w róşnych eksperymentach pomiarowych przy tej samej niepewności uB mogą mieć róşne niepewności uA. Dla menzurandu wejściowego X współczynniki korelacji typu A moşna znaleźć statystycznie [12], ale eksperyment pomiarowy trzeba zorganizować tak, by synchronicznie pozyskiwać kolejne wartości obserwacji wszystkich mierzonych parametrów. Niepewności typu B i związane z nimi skorelowanie szacuje się heurystycznie. Jest to zadanie dosyć zróşnicowane w praktyce, gdyş mierzy się róşne wartości badanych wielkości, na róşnych zakresach jednego miernika lub nawet róşnymi miernikami. Róşne są teş dopuszczalne błędy maksymalne tych przyrządów, a więc i róşne wynikające z nich niepewności typu B [8, 20]. Przy szacowaniu skorelowania typu B w jak największym stopniu trzeba wykorzystywać wszelką pozyskaną wcześniej wiedzę o tego rodzaju pomiarach.
v2 $ ' * ' Z,
*
* W modelu matematycznym zastosowanym w GUM [1] przyjęto załoşenie, şe niepewności składowe uA i uB otrzymuje się dla niezaleşnych statystycznie rozkładów odchyleń wartości obserwacji od estymatora wartości mierzonej. Tak teş jest zwykle w praktyce. Opisana wzorem (1) niepewność mierzonej wartości x pojedynczego menzurandu ma więc taką samą postać jak niepewność sumy dwu wielkości o nieskorelowanych wartościach, z których jedna ma tylko niepewność uA, a druga – tylko niepewność uB. Wartość x moşna więc dowolnie podzielić na dwie części, tj. zamodelować ją jako sumę x = xA + xB o wypadkowej wariancji ux2 = uA2 + uB2 , czyli tak, jak dla splotu dwu nieskorelowanych rozkładów o niepewnościach uA i uB. W pomiarach dwuparametrowych (2D) wartości x1, x2 wielkości mierzonych mona przedstawić jako dwie sumy, których składniki mają niepewności tylko jednego typu, A lub B, tj.: x1 = x1A + x1B
(6a)
x2 = x2A + x2B
(6b)
Rys. 2. Powiązania korelacyjne odchyleń dwuparametrowego menzurandu o niepewnościach typu A i/lub typu B Fig. 2. Correlations between deflections from estimators of the 2D measurand with uncertainties of Type A and/or Type B
Elementy multimenzurandu wejściowego X mogą mieć ogólnie róşne stopnie skorelowania typu A jak i typu B. Wpływa to w tych pomiarach na dokładność wyników pośrednio wyznaczanych wielkości wyjściowych Y. Aby z wyników pomiarów 2D wyznaczyć standardowe niepewności uy1, uy2 i ich współczynnik korelacji ry1,2 korzysta się z wektorowego prawa propagacji niepewności (3). W przypadkach, gdy wnioski z pomiarów mają dotyczyć innych warunków pomiaru niş dla danego eksperymentu, innych niş były przy kalibracji przyrządu, bądź podczas uşytkowania badanego obiektu, to obok bieşących rozrzutów wartości obserwacji o niepewnościach standardowych u1A, u2A i ich współczynnika korelacji rA, trzeba uwzględnić teş przypuszczalne wpływy niepewności u1B, u2B i współczynnika korelacji rB Jako wielkości wejściowe moşna teraz przyjąć [x1A, x2A, x1B, x2B]T, a jako wyjściowe [y1, y2]T oraz naleşy załoşyć, şe nie ma korelacji między odchyleniami niejednakowych typów A i B dla kaşdego z elementów menzurandu wejściowego. Moşe zaś wystąpić korelacja x1A i x2A opisana współczynnikiem korelacji rA oraz między x1B i x2B opisana przez rB. Równanie (3) przyjmuje wówczas następującą postać:
Niepewności wielkości x1 i x2 oznacza się jako ux1 i ux2. Spełniają one równania:
ux21 = u12A + u12B
(7a)
ux22 = u 22A + u 22B
(7b)
=
50
P
O
M
I
A
(8)
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
< =)# ? #
Połowa elementów macierzy kowariancji wielkości wejściowych o rozmiarze 4 × 4 równa się zeru. Wynika to stąd, şe korelacja występuje tylko między odchyleniami tworzącymi niepewności składowe jednego typu A oraz/lub B, zgodnie z załoşeniem przyjętym w GUM [1], iş obie składowe są od siebie niezaleşne statystycznie. W szczególnym przypadku dla bezpośrednich pomiarów składowych wejściowej wielkości wektorowej 2D, tj., gdy Y = X, czyli dla y1 = x1, y2 = x2, wszystkie współczynniki czułości są
jednakowe:
i wyjściowa macierz
kowariancji (8) Uy1,2 ≥ Ux1,2 i wynosi ⎥ u12A + u12B U x 1,2 = ⎢ ⎢ Ď Au1Au 2A + Ď Bu1Bu 2B ⎣
Ď Au1Au 2A + Ď Bu1Bu 2B ⎤ u 22A + u 22B
⎼ ⎼ ⎌
(9)
Tabela 1. Współczynniki korelacji Ď x1,2 wynikĂłw dwu pomiarĂłw o róşnych relacjach niepewnoĹ›ci skĹ‚adowych uA oraz/lub u B oraz róşnych współczynnikach Ď A , Ď B Table 1. Correlations coefficients Ď x1,2 of two measurands for different relations of their uncertainty components uA or uB and different correlation coefficients Ď A , Ď B
Niepewności typu A
Wartości
Współczynnik korelacji
Niepewności typu B
Wartości
Współczynnik korelacji
Współczynnik korelacji ry1,2 elementów menzurandu wyjściowego 2D
Ď Au1Au 2A + Ď Bu1Bu 2B
Lp.
u12A + u12B u1A, u2A
rA
u1B, u2B
u 22A + u 22B
rB
Ď A 1 − k12B 1 − k22B + Ď B k1B k2B
1.
u1A, u2A
rA
u1B, u2B
rB = 0
Ď A 1 − k12B 1 − k22B
2.
u1A, u2A
rA = 1
u1B, u2B
rB
Ď B k1B k2B
rA
u1B = u2B = uB
rB
Ď A 1 − kB2 + Ď B kB2
rA
u1B = u2B = uB
rB = 1
Ď A 1 − k12B 1 − k22B + k1B k2B
rA = 1
u1B, u2B
rB
1 − k12B 1 − k22B + Ď B k1B k2B
rA = 0
u1B, u2B
rB = 0
0
rB
0,5(rA + rB)
3.
4.
5.
6.
u1A = u2A = uA
u1A, u2A
u1A = u2A = uA
u1A, u2A
rA
7.
u1A = u2A = uA = u1B = u2B = uB
rA = 0 rA = 1
u1A = u2A = uA = u1B = u2B = uB
rA = 1
W szczegĂłlnym przypadku, gdy rA = 0, rB = 1 i
(
)
rB = 1 ry1,2 = ½ rB = 0 rB = 1
ry1,2 = 1
otrzymuje siÄ™ ux 1 = ux 2 = 2uB oraz rx1,2 = 1/2
51
W ) ^
) ) & ; ^H U& &&&
W otrzymanej w (9) macierzy Ux1,2 występują niepewności zgodne z przyjętą regułą sumowania wariancji, czyli dla kwadratów niepewności typu A i typu B w obu pomiarach. Współczynnik korelacji między wielkościami x1 i x2 wynosi:
Natomiast warunek
uiB uiA + = kiB + 1 − kiB2 = 1 jest speĹ‚niony uxi uxi
tylko dla kiB = 0 lub kiB = 1, i = 1, 2.
Ď x 1,2 =
Ď Au1Au 2A + Ď Bu1Bu 2B ux 1ux 2
=
Ď Au1Au 2A + Ď Bu1Bu 2B u12A + u12B
Niektóre szczególne przypadki współczynnika korelacji rx1,2 dla wyników eksperymentów pomiarowych o róşnych wariantach par niepewności składowych u1A, u2A oraz u1B, u2B i współczynnikach korelacji rA i rB podano w tabeli 1. Zawiera teş ona skrajne przypadki dla kombinacji wartości rA = (0, 1); rB = (0, 1).
(10)
u 22A + u 22B
Pomiary menzurandu wejściowego X naleşy traktować jako szczególny przypadek wieloparametrowych pomiarów opisanych wzorami (2) i (3), gdy są one powiązane ze sobą tylko w samym obiekcie pomiarowym. Jeśli dla stosunków niepewności składowych i niepewności standardowych ux1, ux2 przyjmie się oznaczenia:
0 ≤ k1B ≥
u1B ≤ 1, ux 1
0 ≤ k 2B ≥
4. Wykresy wypadkowego współczynnika korelacji ry1,2 Ze wzglÄ™du na niezmienność postaci funkcji (11) przy zamianie wielkoĹ›ci u1A na u1B, u2A na u2B oraz rA na rB, wystarczy przeanalizować zaleĹźnoĹ›ci dla jednego typu niepewnoĹ›ci A lub B, gdyĹź dla drugiego typu, przy tym samem zestawie zmiennych, bÄ™dÄ… podobne zaleĹźnoĹ›ci wskutek symetrii tego wzoru. Do przedstawienia współczynnika korelacji skorzysta siÄ™ z ogĂłlnej funkcji wg wzoru (11). Na rysunku 3 podano wykresy 3D wartoĹ›ci wypadkowego współczynnika korelacji wielkoĹ›ci wyjĹ›ciowych jako funkcje Ď y 1,2 = f k12B , k22B i trzech par współczynnikĂłw skorelowania typu A i B dwuelementowego menzurandu wejĹ›ciowego X, tj.: rA = 0, rB = 1; rA = 1, rB = 1 i rA = 1, rB = 0. SÄ… to powierzchnie krzywoliniowe. Zaznaczono na nich Ĺ›lady przekrojĂłw dla trzech wartoĹ›ci współczynnika k2B. Podobne powierzchnie przedstawiono na rys. 4, teĹź dla trzech wartoĹ›ci par współczynnikĂłw rA, rB, w tym dwu ujemnych, tj.: rA = –1,0, rB = 1; rA = –0,5, rB = 1 oraz rA = 0,0, rB = 1. Zaznaczono teĹź jeden przekrĂłj. Wszystkie trzy przekroje z rys. 3 i jeden z rys. 4 przedstawiono nastÄ™pnie na rysunkach 5a–d jako wykresy 2D zaleĹźnoĹ›ci ry1,2
u 2B ≤ 1, ux 2
to otrzyma siÄ™:
u1A ≥ 1 − k12B , ux 1
u 2A ≥ 1 − k22B ux 2
(
i prostszy wzĂłr (10)
Ď x 1,2 = Ď A 1 − k12B 1 − k22B + Ď B k1B k2B
(11)
Dla obu wielkości mierzonych pośrednio współczynnik rx1,2 w (11) zaleşy wyłącznie od stosunków niepewności kiB, gdyş dla 2
2
⎛u ⎞ ⎛u ⎞ kaĹźdego z nich zachodzi ⎜ iB âŽ&#x; + ⎜ iA âŽ&#x; = kiB2 + 1 − kiB2 = 1. âŽ? uxi ⎠âŽ? uxi âŽ
(
)
= 0.25
)
= 0.5
=1,
= 0.81
=0,
=1,
Rys. 3. ZaleĹźnoĹ›ci współczynnika korelacji Ď y1,2 = f(k1B , k2B) menzurandu X jako wykresy 3D dla trzech par współczynnikĂłw skorelowania typu A i B: Ď A = 0, Ď B = 1; Ď A = 1, Ď B = 1 oraz Ď A = 1, Ď B = 0 Fig. 3. Relations of correlation coefficient Ď y1,2 = f(k1B , k2B) of measurand X as 3D charts for three pairs of correlation coefficients type A or B components: Ď A = 0, Ď B = 1; Ď A = 1, Ď B = 1 oraz Ď A = 1, Ď B = 0
52
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
< =)# ? #
Rys. 4. ZaleĹźność w rzucie trĂłjwymiarowym współczynnika korelacji Ď y1,2 dwĂłch wielkoĹ›ci wyjĹ›ciowych w funkcji k1B (lub jego kwadratu) i k2B (lub jego kwadratu) oraz okreĹ›lonych wartoĹ›ci współczynnikĂłw korelacji wielkoĹ›ci wejĹ›ciowych: Ď A = –1,0, Ď B = 1; Ď A = –0,5, Ď B = 1; Ď A = 0,0, Ď B = 1 Fig. 4. Dependences in three dimensional projections of correlations coefficient Ď y1,2 between output quantities as function of k1B and k2B and defined values of correlations coefficients of input quantities Ď A = –1,0, Ď B = 1; Ď A = –0,5, Ď B = 1; Ď A = 0,0, Ď B = 1
a)
k2B = 0,50 (50%),
b)
k2B Âť 0,71 (71%),
= 0,5 (50%)
c)
k2B = 0,9 (90%),
d)
k2B =
1
,
8 0,8
= 0,81 (81%)
= 0,5 (50%),
= 0,25 (25%)
= 1, = 1 = 0,5 , = 1
= 0, = 1
0,6 0,4 0,2 0 -0,2 -0,4
0
10
20
30
40
50
60
70
80
90 100 ,%
-0,6 -0,8 -1
, %
Rys. 5. ZaleĹźnoĹ›ci współczynnika korelacji Ď y1,2 dwĂłch wielkoĹ›ci wyjĹ›ciowych w funkcji k1B (lub jego kwadratu) dla trzech par współczynnikĂłw korelacji rozkĹ‚adĂłw odchyleĹ„ tworzÄ…cych skĹ‚adowe niepewnoĹ›ci typu A oraz B: Ď A, Ď B = {0, 1; 1, 0; 1, 1} oraz róşnych wartoĹ›ci parametrĂłw k2B (i jego kwadratu) Ď A = –1,0, Ď B = 1; Ď A = –0,5, Ď B = 1; Ď A = 0,0, Ď B = 1 Fig. 5. Correlation coefficient Ď y1,2 between output quantities as function of k1B for defined values of correlation coefficients of deflections, which create input uncertainty components type A or/and B: Ď A = 0 or 1 and Ď B = 0 or 1, and various values of k2B
53
W ) ^
) ) & ; ^H U& &&&
⎛u od k12B = ⎜ 1B ⎜u � y1
działania multiplikatywne (mnoşenie, dzielenie itd.), bardzo uşyteczne są wzory wyraşone we względnych niepewnościach. Otrzymuje się je z wzorów dla niepewności bezwzględnych przez podstawienia: ui = xi uri, uAi = xiurAi, uBi = xi urBi dla i = 1, 2. Ze wzoru (10) dla współczynnika korelacji otrzymuje się inną jego postać podaną we wzorze (13):
2
⎞ âŽ&#x;âŽ&#x; , tj.: âŽ
( )
Ď y 1,2 = f k12B
(12a)
Ď y 1,2 = f (k1B )
(12b)
uA2 1 + uB2 1 2 2 = urA 1 + urB 1 x12
(14a)
ur2 (x 2 ) ≥ urx2 2 =
uA2 2 + uB2 2 2 2 = urA 2 + urB 2 x 22
(14b)
Niepewności względne urx1, urx2 jako funkcje niepewności względnych urA1, urA2, urB 1, urB 2 mają postać:
A
R
Y
•
A
U
T
O
M
2 2 urx 2 = urA 2 + urB 2
(15b)
Ď AurA1urA2 + Ď BurB 1urB 2 urx 1urx 2
(16)
Jeśli niepewności względne wielkości wejściowych są znane, np. jednakowe dla całego zakresu, to moşna się bezpośrednio posłuşyć ich wektorowym równaniem propagacji, o strukturze podobnej jak dla niepewności bezwzględnych, podanym w [16], tj.
U δY = Sδ ⋅U δ X ⋅ SδT
(17)
⎥ x ∂y j ⎤ gdzie: UdX, UdY, Sδ = ⎢ i ⎼ – macierze kowariancji dla nie⎣⎢ y j ∂x i ⎌⎼ pewnoĹ›ci wzglÄ™dnych i macierz czuĹ‚oĹ›ci o elementach oznaczonych indeksami i = 1, 2, 3 dla wierszy oraz j = 1, 2, 3 dla kolumn.
Dla niepewności względnych typu A i B moşna w podobny sposób jak poprzednio dla uA i uB stworzyć macierz UdAB, aby otrzymać macierz UdX = SdAB×UdAB × o wyodrębnionych składowych niepewności względnych.
Przy opisie niepewności dla zakresów pomiarowych przyrządów i systemów oraz w pomiarach, których wzory zawierają
I
(15a)
Ď x 1,2 =
X2 $ ' * # 0
M
2 2 urx 1 = urA 1 + urB 1
Współczynnik korelacji rx1,2 w macierzy kowariancji dla niepewności względnych wyraşa wzór zgodny z wzorem (13), a po skorzystaniu z wzorów (15a) i (15b) otrzymuje się wzór (16):
otrzymuje siÄ™, gdy naprzemienne iloczyny niepewnoĹ›ci typu A dla pierwszego menzurandu i dla drugiego menzurandu typu B i odwrotnie, w obu pomiarach sÄ… sobie rĂłwne. Współczynnik korelacji roĹ›nie do 1 dla k1B < k2B, natomiast zmniejsza siÄ™ dla k1B > k2B. Dla wartoĹ›ci k12B < 1 − k22B obserwujemy duşą ujemnÄ… korelacjÄ™ dla krzywych rA = (–1, –0,5), rB = (0, 1).
O
(13)
ur2 (x1 ) ≥ urx2 1 =
u1B u 2B . Zatem maksymalną wartość współczynnika korelacji = u1A u 2A
P
2 2 2 2 urA urA 1 + urB 1 2 + urB 2
W podobny sposób wyznacza się niepewności względne wielkości wyjściowych:
Wnioski szczegółowe wynikajÄ…ce z analizy wykresĂłw na rysunkach 1–5 : − najwiÄ™kszy współczynnik korelacji w caĹ‚ym zakresie k1B ma krzywa rA = 1, rB = 1 gdy k1B = k2B; − dla wartoĹ›ci k12B < 1 − k22B krzywa dla dominuje nad krzywÄ… rA = 0, rB = 1 zbliĹźa siÄ™ do krzywej rA = 1, rB = 1 dla najmniejszych wartoĹ›ci rA = 1, rB = 0 k22B ; − dla wiÄ™kszych wartoĹ›ci k12B > 1 − k22B krzywa rA = 1, rB = 0 jest poniĹźej krzywej rA = 0, rB = 1; − w punkcie k12B = 1 − k22B krzywe dla rA = 1, rB = 0 i rA = 0, rB = 1 przecinajÄ… siÄ™. Współczynnik korelacji wynosi Ď y 1,2 = k1B 1 − k12B i osiÄ…ga w tym punkcie wartość maksymalnÄ… dla k12B = 1 / 2 wartość Ď y 1,2 = 1 / 2 ; − dla wartoĹ›ci współczynnika korelacji A < 0 wystÄ™puje opcja uzyskania ujemnego współczynnika korelacji y1,2 dla fragmentu zakresu zmiennoĹ›ci k1B; − dla rA = –1, rB = 1 i k12B < 1 − k22B otrzymuje siÄ™ ujemny współczynnik korelacji, a dla k12B > 1 − k22B współczynnik korelacji jest dodatni. NajwiÄ™kszy zakres zmiennoĹ›ci k1B z ujemnym współczynnikiem korelacji wystÄ™puje dla maĹ‚ych wartoĹ›ci zakresu k2B. Z dotychczasowych rozwaĹźaĹ„ wynikajÄ… podane poniĹźej wnioski ogĂłlne o szacowaniu elementĂłw macierzy kowariancji w pomiarach poĹ›rednich wieloparametrowych o skorelowanych wielkoĹ›ciach wejĹ›ciowych. Korelacja miÄ™dzy wartoĹ›ciami wynikĂłw pomiarowych dwu skojarzonych ze sobÄ… wielkoĹ›ci, w tym pomiary tej samej wielkoĹ›ci mierzonej dwoma róşnymi miernikami lub na dwĂłch róşnych zakresach tego samego miernika, czyli o róşnych wartoĹ›ciach uB, zaleĹźy od wartoĹ›ci tej niepewnoĹ›ci i jest tym wiÄ™ksza, im wiÄ™ksza jest niepewność pomiaru typu B w stosunku do niepewnoĹ›ci typu A w jednym lub obu wykonanych eksperymentach pomiarowych. Maksymalny współczynnik korelacji ry1,2 = 1 osiÄ…ga siÄ™ dla wielkoĹ›ci w peĹ‚ni skorelowanych skĹ‚adowych typu A i typu B, rA = 1, rB = 1, gdy k1B = k2B. Prowadzi to do warunku
54
Ď AurA1urA2 + Ď BurB 1urB 2
Ď x 1,2 =
Na dole kaşdego z rysunków 5a–d dodano dodatkową nieliniową skalę dla funkcji liniowej
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
< =)# ? #
V2 $ ' * Z,
' ' # %
Wyznaczymy niepewności menzurandu Y = [y1, y2]T otrzymywanego pośrednio z przetwarzania menzurandu X = [x1, x2]T za pomocą dowolnej funkcji Y = F(X). Wartości i macierz kowariancji menzurandu X są znane, np. otrzymano je z bezpośrednich pomiarów. Wartości menzurandu wyjściowego Y = [y1, y2]T będą miały niepewności uy1, uy2. Oba menzurandy są dwuparametrowe (2D). Do wyznaczenia niepewności skorzysta się z ogólnego wektorowego równania propagacji niepewności (3). Funkcję przetwarzania linearyzuje się do postaci macierzy czułości S funkcjonału F( ) o wymiarach [2 × 2] i otrzymuje się:
obserwacji pomiarowych moşe być wymagane wyodrębnienie niepewności składowych typu A i B wielkości menzurandu wyjściowego Y oraz skorelowania dla par tych wielkości. Zanalizujemy to na przykładzie menzurandów 2D w oparciu o zaleşności (7a) i (7b) dla X = [x1, x2]T, tj.: ux21 = u12A + u12B i ux21 = u22A + u22B . Wejściowa macierz kowariancji UX jest sumą dwu macierzy kowariancji dla składowych niepewności typu A i typu B wielkości wejściowych X: UX = UXA + UXB
(22)
dla niepewnoĹ›ci typu A ⎥ u12A U XA = ⎢ ⎢ Ď Au1Au 2A ⎣
(18)
⎥ u12B U XB = ⎢ ⎢ Ď Bu1Bu 2B ⎣
∂y1 ⎤ ⎼ ∂x 2 ⎼ ⎥ι1 Îą 2 ⎤ ⎢ ⎼; ⎼≥ ∂y 2 ⎼ ⎢⎣ β1 β2 ⎼⎌ ∂x 2 ⎼⎌
u 22A
⎼ ⎼ ⎌
(22a)
oraz dla niepewności typu B
gdzie: macierz funkcji F()
⎥ ∂y1 ⎢ ⎢ ∂x1 S =⎢ ⎢ ∂y 2 ⎢⎣ ∂x1
Ď Au1Au 2A ⎤
Ď Bu1Bu 2B ⎤ u 22B
⎼ ⎼ ⎌
(22b)
Transformację macierzy kowariancji UXA i UXB składowych niepewności menzurandu X realizuje się po linearyzacji funkcjonału Y = F(X) wg wzoru (2) w następujący sposób:
a1, a2, b1, b2 − odpowiednie współczynniki.
UYA = SU AS T
(23a)
Po wykonaniu działań wg wzoru (18) otrzymuje się niepewności wyjściowe
UYB = SU BS T
(23b)
uy21 = Îą12ux21 + Îą 22ux22 + 2 Ď x 1,2Îą1Îą 2ux 1ux 2
(19a)
uy22 = β12ux21 + β 22ux22 + 2 Ď x 1,2 β1β 2ux 1ux 2
(19b)
gdzie S dane jest zaleşnością (3a).
oraz współczynnik korelacji ry1,2 dla zmiennych y1, y2:
Ď y 1,2 =
Îą1β1ux21 + Îą 2 β2ux22 + (Îą1β2 + Îą 2 β1 )Ď x 1,2ux 1ux 2 ux 1ux 2
Macierz kowariancji menzurandu wyjściowego UY moşna wyznaczać dwojako: oszacować macierz UX i z niej wyznaczyć bezpośrednio UY, lub z macierzy UXA i UXB znaleźć macierze UYA i UYB dla obu składowych niepewności uyi wszystkich wielkości wyjściowych. Otrzymuje się: UY = SUXST = S(UA + UB)ST
(20)
(24)
gdzie: Z (20) wynika, şe rx1,2ux1ux2 = rAu1Au2A + rBu1Bu2B. Wykorzystuje się to w dalszej analizie. Z (7a, b) i (19 a, b) wynikają równania dla niepewności wyjściowych uy1 i uy2.
⎥ u12A + u12B UX = ⎢ ⎢ Ď Au1Au 2A + Ď Bu1Bu 2B ⎣
uy21 = Îą12 (u12A + u12B ) + Îą 22 (u 22A + u 22B ) + 2Îą1Îą 2 ( Ď Au1Au 2A + Ď Bu1Bu 2B )
Ď Au1Au 2A + Ď Bu1Bu 2B ⎤ u 22A + u 22B
⎼ ⎼ ⎌
oraz (21a) UY = UYA + UYB = SU AS T + SU BS T uy22 = β12 (u12A + u12B ) + β 22 (u 22A + u 22B ) + 2 β1β 2 ( Ď Au1Au 2A + Ď Bu1Bu 2B )
(21b)
_2
* (
' * Przy stosowaniu przyrządów i układów pomiarowych w róşnych warunkach otoczenia i przy róşnym losowym rozrzucie
(25)
Jeśli wyznacza się UY, to oba sposoby są równowaşne. Sprawdziliśmy to dla modelu o dwuwymiarowych (2D) menzurandach X i Y. Ale jedynie z wzoru (25) niespotkanego w literaturze wyznaczyć moşna niepewności składowe elementów menzurandu wyjściowego UY i współczynniki korelacji dla kaşdego z ich typów A i B. Ogólnie dla splotów rozkładów wielowymiarowych tworzy się sumę macierzy kowariancji przy oddzielnej transformacji ich niepewności typów A i B.
55
W ) ^
) ) & ; ^H U& &&&
‚2 & 0* Y2
8.
W części 1 przedstawiono propozycję wektorowej metody wyznaczania niepewności w pomiarach pośrednich wieloparametrowych. Jest ona rozszerzona w stosunku do podanej w Suplemencie 2 do Przewodnika wyraşania niepewności GUM [1, 3] i uwzględnia przypadek skorelowania wielkościach wejściowych. Podano wzory dla niepewności i współczynnika korelacji na wyjściu przy znanych niepewnościach składowych A oraz B wielkości mierzonych i znanych współczynnikach skorelowania kaşdej z nich z jej odpowiednikami dla pozostałych wielkości. W zastosowanym tu modelu załoşono, şe dla kaşdej z wielkości mierzonych jej składowe niepewności typu A oraz B są od siebie niezaleşne statystycznie i sumują się geometrycznie (pierwiastek z sumy ich kwadratów). Stąd wynika, şe w pośrednich pomiarach wieloparametrowych skorelowanie mierzonych wielkości wejściowych jest tylko rezultatem osobnego skorelowania odchyleń od estymatora wartości tworzących niepewności A oraz B. Umoşliwia to oszacowanie wpływu tych wszystkich skorelowa na wartości elementów macierzy kowariancji. Macierz ta jest podstaw oceny dokładności przy łącznym dalszym stosowaniu wielkości ze sobą skojarzonych. Współczynnik skorelowania między bieşącymi odchyleniami tworzącymi niepewności typu A dwu mierzonych wielkości moşna wyznaczyć eksperymentalnie, jeśli przeprowadzi się ich pomiary synchronicznie. Natomiast współczynnik skorelowania między zbiorami odchyleń opisanych niepewnościami typu B ocenia się heurystycznie w oparciu o dotychczasową wiedzę opartą na doświadczeniu z wielu innych eksperymentów i wnikliw analiz danego eksperymentu. Zastosowanie w praktyce przedstawionej rozszerzonej wektorowej metody propagacji niepewności przyblişy omówienie kilku charakterystycznych przykładów. Będą one analizowane w części 2. i dotyczą estymacji niepewności oraz wypadkowego współczynnika korelacji dla przetwarzania wyników w pośrednich pomiarach dwuwymiarowych (2D) według funkcji liniowej oraz kilku rodzajów funkcji nieliniowych. Wyniki osiągnięte w tych przykładach uzasadniają, by zaproponowaną tu metodę stosować w praktyce pomiarowej, w tym przemysłowej. Zostaną teş podane łącznie wnioski końcowe, wynikające z obu części pracy.
9.
10.
11.
12.
13.
14.
15.
16.
17.
~ # % 1. BIPM, IEC, IFCC, ILAC, ISO, IUPAC, IUPAP and OIML, Evaluation of measurement data — Guide to the Expression of Uncertainty in Measurement, JCGM 100:2008, GUM 1995 with minor corrections. 2. Supplement 1 to the ‘Guide to the Expression of Uncertainty in Measurement’ – Propagation of distributions using a Monte Carlo method. JCGM 101:2008, BIPM. 3. Supplement 2 to the ‘Guide to the Expression of Uncertainty in Measurement’ – Extension to any number of output quantities. JCGM 102:2011 BIPM. 4. Supplement 4 to the ‘Guide to the Expression of Uncertainty in Measurement’ – An introduction to the “Guide to the expression of uncertainty in measurement� and related documents JCGM 104:2009, BIPM. 5. Supplement 6 to the ‘Guide to the Expression of Uncertainty in Measurement’ The role of measurement uncertainty in conformity assessment JCGM 106:2012, BIPM. 6. Conformity Assessment – General Requirements for the Competence of Testing and Calibration Laboratories, ISO/ IEC 17025:2017. 7. EA-4/02.M: 2013 Evaluation of the Uncertainty of Measurement in Calibration.
56
P
O
M
I
A
R
Y
•
A
U
T
O
M
18.
19.
20.
21.
A
T
Warsza Z.L., Metody rozszerzenia analizy niepewnoĹ›ci pomiarĂłw. Monografia, Oficyna Wydawnicza PIAP, Warszawa 2016, ISBN 978-83-61278-31-3. Zakharow I.P., Estimating measurement uncertainty on the basis of observed and logical correlation. “Measurement Techniquesâ€?, Vol. 50, No. 8, 2007, 808–816. Dorozhovets M., Warsza Z.L., Udoskonalenie metod wyznaczania niepewnoĹ›ci wynikĂłw pomiaru w praktyce. „PrzeglÄ…d Elektrotechnicznyâ€?, R. 83, Nr 1, 2007, 1–13. Dorozhovets M., Warsza Z.L., Propozycje rozszerzenia metod wyznaczania niepewnoĹ›ci wyniku pomiarĂłw wg Przewodnika GUM (2) UĹ›ciĹ›lenie metod obliczeĹ„ niepewnoĹ›ci typu B. „Pomiary Automatyka Robotykaâ€?, Nr 2, 2007, 45–52. Warsza Z.L., ZiÄ™ba A., Niepewność typu A pomiaru o obserwacjach samoskorelowanych. „Pomiary Automatyka Kontrolaâ€?, R. 58, Nr 2, 2012, 157–161. Warsza Z.L., Ezhela V.V., Zarys podstaw teoretycznych wyznaczania i numerycznej prezentacji wynikĂłw pomiarĂłw poĹ›rednich wieloparametrowych. „Pomiary Automatyka Kontrolaâ€?, Vol. 57, Nr 2, 2011, 175–179. Warsza Z.L., Ezhela V.V., O wyraĹźaniu i publikowaniu danych pomiarĂłw wieloparametrowych – stan aktualny a potrzeby. „Pomiary Automatyka Robotykaâ€?, Nr 10, 2011, 68–76. Warsza Z.L., Evaluation and Numerical Presentation of the Results of Indirect Multivariate Measurements. Outline of Some Problems to be Solved. Advanced Mathematical & Computational Tools in Metrology and Testing IX, ed. by Franco Pavese, Markus Bär et all, Series on Advances in Mathematics for Applied Sciences, Vol. 84, World Scientific Books 2012, Singapore, 418–425. Warsza Z.L. Part 1, and Warsza Z.L., Puchalski J. Part 2: Estimation of uncertainty of indirect measurement in multi-parametric systems with few examples. ppt: in CD Proceedings of conference: Problems and Progress of Metrology ppm’18 Szczyrk 04-06. June 2018. Series: Conferences No. 22, Metrology Commission of Katowice Branch of the Polish Academy of Science. Warsza Z.L., Puchalski J., Estimation of vector uncertainties of multivariable indirect instrumental measurement systems on the star circuit example. XXII World Congress IMEKO 2018 Belfast. Abstract in CD Proceedings PO-062 and IOP Conf. Series: Journal of Physics: Conf. Series 1065 (2018) 052026, DOI: 10.1088/1742-6596/1065/5/052026. Warsza Z.L., Puchalski J., Estymacja macierzowa niepewnoĹ›ci wieloparametrowych pomiarĂłw poĹ›rednich z przykĹ‚adami. „Pomiary Automatyka Robotykaâ€?, R. 22, Nr 2, 2018, 31–39 DOI: 1014311/PAR_228/31. Warsza Z.L., Puchalski J., Ocena niepewnoĹ›ci prostokÄ…tnych skĹ‚adowych impedancji wyznaczanych poĹ›rednio z pomiarĂłw skĹ‚adowych biegunowych i vice versa. „Pomiary Automatyka Robotykaâ€?, R. 22, Nr 3, 2018, 61–67, DOI: 1014311/PAR_229/61. Warsza Z.L., Puchalski J., Estymacja niepewnoĹ›ci wartoĹ›ci funkcji z pomiarĂłw w punktach kontrolnych. „Pomiary Automatyka Robotykaâ€?, R. 22, Nr 4, 2018, 39–50, DOI: 1014311/PAR_230/39. ZiÄ™ba A., Analiza danych w naukach Ĺ›cisĹ‚ych i technice. Wydawnictwo Naukowe PWN, Warszawa 2013.
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
< =)# ? #
( % % B % # ? # B B 0 * . A # % B % V F Abstract: The two-part work presents an extended version of the vector method of uncertainty evaluation of multiparameter measurements given in Supplement 2 to the Guide on Expression of Uncertainty of Measurements about the English acronym GUM [1]. The discussion was made on the example of measurements of a two-parameter (2D) measurand with correlated parameters. A novelty in the proposed method is to consider the correlation of deflections from estimators of values of quantities measured directly on the input, which created individual components of the type A and/or type B uncertainties. In this part 1, general formulas for the covariance matrix, output uncertainties and correlation coefficient were determined. Formulas for several characteristic specific cases are presented in table 1. The graphs show the correlation coefficients of the output quantities as a function of the type B contribution in the full uncertainty of the input quantities. In Part 2 will be presented examples of estimation of uncertainty and resultant correlation coefficient for several types of linear and non-linear functions of processing measured results. The general conclusions will be given on the end. Keywords` & % )6& % b # 8 % A& ( # & b
" " @ = ! B
" ; !
*+ J'(
'( (
8F > % ? ( 0 > *+E+& % = *+JQ& % % *+Q! 0 [ 4 ? *+EKN*+J *++-N*++E& 0 > *+J!N*+Q!& 0 a *+Q!N*+QK L ( % > % . G % M& $ ( [ $ % 8 = 0 4 B = ( ; % > % *+QKN*+K) % 8 0 < 4 0 *+K N*++) 6 % B ?% X % *++)N*++E& 0 G % *+K N)!!) $F ( < 0 4 = 8 0 < 0480 8 % -! F & J ( # & % F % & ** < = ) % < 0 0 ( . B ( ( 0 08G 8 % B (
8F > % / . = B L*+KJ M > % ? = L*+KK M 0 > = > *+KQN*++E 0 > = & *++E F % S " *! F @ ( # = % S S= S 0 ( = # % B / < % S% T # $% )!!J ; < % B % S ( % < < % % = ( & F ( < ( % =
57
NR 3/2015
58
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
) # ) " $
K((% UPOT cUOF "& Oh % UBOGUc
Informacje dla Autorรณw < # # A =
)# ) % #
( ? ! = Oc )
OGUh & M $#
) # A ! #) H
#
# ] =
$ ' ) )
) ) #$
) ] #
#$
) # M
E ]
] $
) $ '
$ =
# $ = =
)
# Z$ $ = E [& * # M
# ! #$
) #
) = E )
) ) #$
= ] )
M ! = #$ )
&
Wskazรณwki dla Autorรณw ' # (
' # # =
#$ # # ) Pomiary Automatyka Robotyka
H # ]
) + A # # # Z ]
QG M [ #
)
= ) A ) # B # M
) E Z #
#$
#
[ A # # Z $ ^
Uย GAOGG M [ #
)
= ) A # Zย AQ [ #
) #
= ) = ) A ^H # # A #
) Z #$ & = )[ A
# ) #
) # = ) A # $ #
# = ) A # B= E B $
) & & & = #$ & a
^ ) & hGG ) &
UGGG
^
= #) ] = E ^ &
#
) H $ ^H M ]
) G F # =
H $ ^
U # = ZPG GGG M
) #$ hGGG )O #
M [ & Q ! =
#& ?
# # M
] $ ^H #=
^ & Nie drukujemy komunikatรณw! * # # )# ] $ H ]
$ = E # M
Z # # #[ $ ^
ย GGATย G M # = E & " = $
M M
)
) = &
#$ ] # # ]
' ]
H # #&
&.! ! ' (
( # ย *
'
' ย ' ( '
&
#
] ]
^H ^H # #&
# ) M
#$ ]
=
# = &
< M ^H
# $ )# ] H
$ #
)
)# ] $ H =
)
)
l M &
7 * ' .
! !
M = # &
M ^ & ]
# ! # # ]H
) # ]
# # A # M E
$ = E & *
#=
" M
# ) =
) ) & ?
#
# !
$ )M
)
" ]&
Kwartalnik naukowotechniczny Pomiary Automatyka Robotyka jest indeksowany w bazach BAZTECH, Google Scholar oraz INDEX COPERNICUS L4 V )!*J[ Q &)-M& @ w bazie naukowych F @ ARIANTA. Punktacja MNiSW F K L *))-M 0 S S % % $ X & % F F w kwartalniku naukowotechnicznym Pomiary Automatyka Robotyka.
ย c
K%i " ; D * WC "Â??
7*
* %
'
' 0
' ' " # = Pomiary Automatyka Robotyka
! ]
% #
( ? ! = ] ]
#) =
# M ) =
# B # M # M
)
# # # # # #
#
^ ] = + 1. % 0
U #
wymieniowego Autora ' # ' A
) #
M
& _=# # ` &
# M
#$ M$ M
#
#
#$ $
) $ =M
) ) & 2. % 0 # 0
' Â
( '(
jej powstanie A ) + A !
$ ) ] #
#
# = $
M # $ $
M ) # B # # ]
)
#
#$ A ! #
_= =` &
# # M
$ ]
# #
) M #
)
)
# $
M =
#
#$ Â?
FG
P
O
M
I
A
" ) #
)
# #
] ) #
] ] ) !
] )
= )
) $# ) &
3. % 0 ) (
<Â
'
' A
#
E
#$
# # $
#$
)
) =
)
l M ^ M
! &
_E # ` A
) $ =
# )
) ) ' #
#$
) #$
E $ '&
/ ' ' (' # ) # )
#$ # ) M ! )
)# ] #
$
$ #
# !
#
A ) ) ) !
H = H ^
# & < ) B
^
]! # &
„
'(
przeniesienie praw
:: ~Â… < ]
# #
#$ #
) ]
#
& W)
) # #
& ?
# =
# # # #)
# = ] #&
Redakcja kwartalnika Pomiary Automatyka Robotyka '
/ # M
#$ # # )
) # ) " $
)
# M
] ) #) * K &+
" # $ % & B 0 % % & . / $ % 046 # / % 8 9 ( & Œ0 Automatyka Robotyka�, ISSN 1427-9126, R. 23, Nr 1/2019, 5–10, DOI: 10.14313/PAR_231/5.
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
czasopisma
pomiary
sprawdzian
miara
POLSPAR
eksperyment
automatyka PIAP
#
seminarium
kalendarium
szkolenie
kwartalnik
federacja
nauka
publikacje
automatyka
stowarzyszenie
HORIZON 2020 ]! innowacje organizacja projekt konkurs
konferencje
relacja
POLSPAR
POLSA
publikacje
AutoCAD streszczenie
agencja kosmiczna
dr h.c.
innowacje
IFAC
ZPSA
l
profesura
recenzje
relacja
szkolenie
doktorat
robotyka seminarium
sterowanie
S
esa
szkolenie
62
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR
1/2019
D; %D (K’“ K
B ( O " 8 G F 0 8 % X .
63
NASZE WYDAWNICTWA
www.piap.pl
64
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
b >K KC ; j ? %( % W ?˜
! " " # $ %
> % ( ( $ ) 6 $ @ @$ $ ! " 4 $ $ #
0
) $ ! 0 $ ! 0 $ 0* # & #
& ' ! " .
#
$ = ]=
! ) = E Wielooperacyjne ste " ' ! & ' )
]
)
#$ # M
) #$ ) E E ) # # Z $ ? ([& @ M ) )
# $ $ M ) # ) ) ) ) )
^
# ) ] ) = )M
# ) = # )
$ = ) & $
= ] = M ]
) )
)
# ) )
# ] # )
# ^ ]= ) # = & ;
C ) ? # )
# # # #
] ? ) @ OGUQ K(>% cTQ Qh TPOU OQÂ P & OGG&
$ '
] # ] '+ U&
)
#
# ) # M! $ ! = ' # ] M
$ =
#&
O&
= # #
# ) # # #
] # )
#& h&
( = )M =# M
]
$M
# # #& P&
* $M # $ =
) M ) )
^ # # # E
)
M
! # # #& W # ) $
E
^ = = # #
# ) # ) = # ) =
= ] # $ ' ) ) ) # E =
#& >
) # E ) $ > ( ) ^ #
( ) ( ) )
K & ) ) E $ + U&
) ) ) # M
$ ' )#
) ! ) )
) # ) ) ) ) #& O&
$ ' )
b
]
^ = ) #
= =
$ ' ) & h&
? =M
$ ' ) ] ^H = = # #
) = ]
! " " # $ % R F > % ? ( 0 ;% T > )!!* & > % $ $ F % S % # "
" ronowych& < F % F P B )J )!*+ T % F ( F & I & % G % X S > % ? 8 0 ;% T $% T % < *++! % F " % 8 $ > @ B & 8 % B & F B ( ;% % S ( % T S & (< S S F % ( < % F ( % S
# # = M $ $ M!
)$ ) ! ) $M # # #
] ' =# M &
B " T
FÂ
KALENDARIUM
Kalendarium wybranych imprez Nazwa konferencji
Data konferencji C ) = '
Informacje dodatkowe
15th IFToMM World ; = OGUc
hG&GFAP> B OGUc
Uย B GU B OGUc
Krakรณw
www: +BB ย ))OGUc& )
mail: = V ย ))OGUc& )
hAย B GT B OGUc Uc B UUB OGUQ
@ '
www: +BB & # =&= & B XOGUc
mail: XOGUcV # =&= &
UOth K ?
on Robot Motion and Control " ; OGUc
GQAUG B GT B OGUc
'
www: +BB ) & # & &
mail:
&) V # & &
5th IFAC Conference on Intelligent Control and # ) ( K; %( OGUc
OUAOh B GQ B OGUc P B GO B OGUc
>
Wielka Brytania
Qth Ki ; ( ) #)
( )
D;b C" %K;( OGUc
GPAGF B Gc B OGUc
U B GF B OGUc
? ' #
www: +BB &) OGUc& = mail: V) OGUc& =
Tth Ki ; ( ) #) ( )
( # # ;
(((; OGUc
GcAUU B Gc B OGUc
U B Gh B OGUc
Sinaia "#)#
www: +BB & &# X& B OGUc
mail:
X X V # ) &# X&
UPth Ki ; ( ) #)
* =
DX # b#)
( )
b ( OGUc
UFAUc B Gc B OGUc
UF B Gย B OGUc
Tallinn D
www: +BB & #& B ) OGUc
mail: Ki ;b (OGUcV #&
K #
Konferencja Metrologรณw
OhAOย B Gc B OGUc
Opole
www: +BB) )OGUc& & & B
mail:
)& V & &
UFth IFAC Conference on = )) $ * X
D)$ ( )
* (OGUc
OcAhU BUG BOGUc
C ? e (
15th D# ?
X ;
* =
;* OGUc
OUAOO B UU B OGUc
hU B Gย BOGUc
Bolonia ?
www:
+BB X &# $ & B OGUc
mail: OGUcV# $ &
OU Ki ; ? ; =
UOAUT B GT BOGOG
Berlin Niemcy
www: +BB & OGOG& =
hG&GQAGh B Gc B OGOU
Jokohama Japonia
UGth Ki ; ( ) #)
K = # ) #
g
K g OGUc
ffKKK ? ; =
K # )
Confederation K D OGOU FF
P
O
M
I
A
R
Y
โ ข
A
U
T
O
M
A
T
Y
www: +BB &ย #$& &# B OGUc mail:
&) Vย #$& &#
www: +BB & #
mail:
& X V #$ &
www: +BB & ) OGOU& =
mail:
) OGOUV & & &
K
A
โ ข
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
%iD"D%; D j
< ?KD*ÂŁ
E* B % " # B (< B ! +E +& $ F +,'- " ( ( ) 3 0 (G 0 $ ( ) *$"
Komitet Organizacyjny MKM 2019 ] ) # % # =
A
& $& !& ? ) ] ) #
= =
A $& !& # " ]
& ? ] ) #
% # =
A $& !& ( ) @ ^
& ;
Konferencja organizowana jest & ! & = A " A
& & $& !& C#
A A
# ? ^ = A ) = %
A ( D M A @ M = W #
# =M
] ] = ] ^
#
# &
] %
D = D = $
^ ( UcFF &
$ ) # ) & ?M
) '
) =M A # %
=M A
^
)
)
) # M
! #!] # ) )
& #= % = & D )# " )
UcFT & @ A
&
UcFQ & #
# )# % )
# )
#
) # )
M # A M$ M = #
H #
# = ) =M & @
M ) # &
( M M
!
) # # &
 G = ^
# # ]
M$ )# ] ) = ]& "M ! Â U&
$ ]
M ^ #
)# ] ) = ] !
) # # )# ] ) ) & =
] #
) = ! M$
) # M (
D M
M
¢ ] ) =
^ ]& $ ) # ^ ) =
] # ]
) + U& % # ) = O& = ^
h& ) ) ) # )
) P& )
) &
> # & -2!
<
0 0 -& # ! =
& ! ,?, @B4E
A l ^H
# ) = Â?
% " 1 - "
A ) = = Â?
% " >& 0
A ) E Â?
K #=#
A ) ^ E ) Â?
#
= $ ]
A ) $ = ) Â?
$ &
A )
) = Â?
$ $ ]
A = M
) Â?
_< ) # $ ` A &
A M
) # ^ Â?
$ $ ]
_< ) #
A ) )M
) Â?
$ ` A &
A ) = &
FT
%iD"D%; D j < ?KD*ÂŁ
XX Krajowa Konferencja Automatyki KKA 2020 ) ) ) $ ) )
( ( $ H $ ! ++ +O $ +,+, " % ) $ 0 1
!
H*
)3 ) 0 )( 0
+, $ +,'U "
Krajowa Konferencja Automatyki jest cyklicznym (organizowanym w cyklu trzyletnim) wydarzeniem naukowym. Konfe = ! ) )
) ^ # # ) )# ] )
= ) # ) &
Stanowi tradycyjne krajowe forum prezentacji wynikĂłw oryginalnych prac badawczych i aplikacyjnych w dziedzinie automatyki, techniki systemĂłw i robotyki.
projektowanie i dobĂłr =# M # M
# = )
# M *
' )
$ # ) & O&O
W
# M
! ) )
) # #
pasywne, liniowe regulatory # & 2.3 Projektowanie sterowania # ) )
parametrach, sterowanie adaptacyjne, sterowanie odporne, ^ = #
# ]
predykcyjne, sterowanie $ ^ # ' M
$ '
# & 2.4 Sterowanie optymalne = M! & O&Â ( # ]
Z X
control), systemy ekspertowe w sterowaniu.
Tematyka konferencji obejmuje 0
` Y2
- Â %
( * *
( ‰ U&U
# )
) ]= ) )& U&O
% # )
) ]= ) )
# & U&h ( ) '
(discrete event systems). U&P ( ) !
# # ) & 1.5 Systemy stochastyczne. U&F W '
) #
) ! & 1.7 Systemy rozmyte i sztuczne sieci neuronowe. U&Q W * #
# = #&
v2
- Â
'
' %
‰ 3.1 Sterowanie komputerowe, mikrokomputerowe, procesory = =#
$# & W #
Z2
- '
( ‰ O&U
W #
! ) )
68
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
Pomiary Automatyka Robotyka, ISSN 1427-9126, R. 23, Nr 1/2019
3.2 h&h
h&P
h&Â
rzeczywistego. Interfejsy ) # &
# # & Systemy cyberfizyczne (cyberphysical systems).
) = # #
inteligencji i inteligencji obliczeniowej, algorytmy heurystyczne, badania ! & W )
H !
# # ) & ( # ]
wykonawcze automatyki.
4. Robotyka ( ‰ 4.1 Konstrukcja, projektowanie i sterowanie manipulatorami i robotami stacjonarnymi, $ ) ) $ ) ] )
] ) & 4.2 Roboty inteligentne, czujniki #
)
wizyjne, fuzja informacji, roboty autonomiczne. P&h K $
$ $
rehabilitacyjne i ratunkowe. P&P ( ) ] $ M
M )#
systemy agentowe. X2
/
' ‰ 5.1 Internet Rzeczy (Internet of C = [ ) P&G
(Industry 4.0) i inne systemy i metodologie sterowania systemami wytwarzania. Â &O
i sterowanie procesami
$ )
$
(e-enterpise, enterprise control systems), systemy sterowania w logistyce. 5.3 Integracja systemĂłw sterowania
$ )
!
komunikacja i systemy sieciowe. V2
' ( ‰ 6.1 Sterowanie procesami ) )
) ^
! ) )
^ M & F&O
(
) ^
wydobywczym, hutnictwie i przetwarzaniu kopalin. 6.3 Systemy sterowania w energetyce, przetwarzanie energii.
_2
' ‰ T&U ( # ) $ ^
M ' # )
# = & 7.2 Systemy sterowania w pojazdach ] )
i diagnostyka, dynamika pojazdu, integracja systemĂłw sterowania
l & > '
i diagnostyka, pojazdy autonomiczne. T&h ) = '
#& T&P W $
]
i podwodnych. T&Â W $
]
i wirnikowych. ‚2
'
# Â
# ‰ Q&U
w rolnictwie. Q&O
w systemach biologicznych, medycznych, w fizjologii, sporcie i psychologii. Q&h
w ekosystemach.
)$ ] = M =
# ) ^ M $
) # # M
blematyki automatyki i robotyki. ? )] #
# !
)
)# ] # ) ] ! )
projektowaniem systemĂłw sterowania. Istotnym efektem, jaki jest planowany do ]=
= M $ ! # # ) # ) ! ^H
) ^ '
= ] M ] M &
#
= szonych referatĂłw w serii Advances in Intelligent Systems and Computing wydawnictwa Springer. ? ) https://www.kka.p.lodz.pl/
Za Komitet Organizacyjny KKA 2020 # % F @ 8 % A & % 08X
\2
' ( ' (
( ‰
c&U
w systemach finansowych, biznesowych, makro i mikroekonomicznych, w systemach
& 9.2 Inteligentne miasta (smart cities). c&h
(
automatyki i robotyki. 9.4 Edukacja w obszarze automatyki i robotyki.
, 0 * ' ' % ( 0
0'
* Warszawa
1958
>
UccU
?
UcFU
@ '
Ucch
Gliwice
1964
Opole
1999
KrakĂłw
1967
Zielona GĂłra
2002
@ '
UcTU
Warszawa
2005
'
UcTP
Szczyrk
2008
RzeszĂłw
1977
Kielce
2011
Szczecin
1980
?
OGUP
JM l
UcQÂ
KrakĂłw
2017
Lublin
1988
69
D; %D (K’“ K
. O G [ B % ( % ? 8 0 $ # ) 0 (
$ # *$ Ethical Aspects of Research and Engineering Methodological and Ethical Aspects of Research3 * +,,V " 0 $ ! B $
)"
X0 0Y 0 $ )
3 ! 0 *$ > [
&
! ' ) badawczej w naukach empirycznych\3 # ( # B $ $ ] $ B $
$ +,'' "
]!
) M A
! ) l ) #= ) A
H ) =] ! # &
# ) )
M ] # ]
! ) # A ) ) =
=M ) = # !
] )
# ] # M # ] #
& ( ) ) )
] ] $ # ) = ] ]
M
) = $ ' #
&
STEM
Roman Z. Morawski
TECHNOSCIENTIFIC RESEARCH METHODOLOGICAL AND ETHICAL ASPECTS
# ! #
^
) =
]
$ )
$ H ] +
) ) & ? M !
# ! H ! $ = # =
# $
#) ) = $ ' # $ '
# & > ! ] ^ H = $ = # ]
Z ) [
)# ]
$ ) Z
) = ] [&
? ] ^H $ =
Œ" S@ # # & @ % Z = # ( % & @ S % % % ( F ( T ( % < ˜ G B & < " J "= Methodological and Ethical Aspects& > % ; & A =A )!*+& [\\ % ( \ \ % \-+K+)J™ A L06/M N =4 AX +QK= =**=!EK-!J=J L = % M& A L?0 AM N =4 AX +QK= =**= !EK-*)=Q L = % M& F N 4 AX +QK= =**= !EK +!=K& E!)
70
P
O
M
I
A
R
Y
•
A
C # ]! ] #
# M # # # # = ) # )
)
#
& %
) ! $
# H ) ) # $ ) ^ + = )
# ) A )
M ) $ '
] $ ] ^H
) #
#
) ) $ #)
$ ' ) & < ! ^H
$ ) #
$ U
T
O
M
A
T
Y
K
A
•
R
O
= )M $ ! M ^ ) =
& > $
^ $ $ $ ) # M
= $ M ) ) $ = ) $ $
$ #) Z
^ ] $ ¢ # ] $ $
# [&
]!
) ^ # ) ) =
M #
^ # = ) ]=#
M =
B
O
T
Y
K
A
NR 1 / 20 1 9
) # ) " $
K((% UPOT cUOF "& Oh % UBOGUc
0 " ! " " @ # $ $
G B S > % ? . 4 # 0 > % *+JQ 0 % % ( % L % 9 M ( *+JK & % F % L% S ( % # ( M <^ ( I% % % F % S ( S[ ( < S F < % <@ F @ & # & & F ( % > ( # F B % % ( < L*+K+ M % % % & < I % % 0 & " % & 9 % > A R F < % !! <
$ ) =
' ! ) # ) $ = ) ! A ^
# = = ! ^ &
) )
$ $ M M! ] ) ] M ]!
) ) = $ ' # M # ) ] # & ; ] =#
= H
$ ' Z [ # H # ย ) H !
]! # ) =] H
$ #
) # ^
) $
) $ ) > =
)
)
^ ยฅ * ^
# ) )
M !
# ! $ ) $
^
= ) #& ]! C " ยฅ ] $ + # # ) = $ '
# # ] ) ) ) = Z ) # = [ ) #
M
^ ) = ! ] ) ) ) )
! =
) ) & ]! )
) !
)
# ย # =# #) ! ] )
)# #
= $#
) = '
= )
) M &
M$]
#
) # )
] $ ] ] # ]+
] = $ $# ]
) & !
M$]
! )
# M !
)
! ] M )
M +
ยฆ ) ]
= ) $ ) A ]! ) ) = # # =M ย
ยฆ ]= M
) ! $ $ $ #
#$ # ! ' $ $ M &
& N F Q! < F 4??? . 4 % B & B & A % ( 0 ( % R # R @ < < S *- < ( L 8& " % & ?M G B % ( %< & ( # 6 N % N % < % & % ( N % ( F % T R F < % < $ F % +! < % < S@ F @ ( R @ & % S > % [ ( & # & & ย # # S S F % T ย (
* = ) ) ) ) Hยง ) =
! )# ) = H = $ ' #
$ ! ^ ) = & ]!
= ] ' #
$ )M =
] & # =
= ]! ! ^
] # ) !
^ H #! + = =M )
) = # ]
# ^
# &
G % % # 8
* = ^ ^ ) # # = ] = ] #
$ )M ) ) ! # # # ]
!
# ) M ) = = ) ! ) H
^ ! =#)
) ] ) ] )
)
A
&
71
NASZE WYDAWNICTWA
VOLUME 13 N°2 2018 www.jamris.org pISSN 1897-8649 (PRINT) / eISSN 2080-2145 (ONLINE)
Indexed in SCOPUS
www.jamris.org
PAR P O M I A RY • A U T O M AT Y K A • R O B O T Y K A
PAR P O M I A RY • A U T O M AT Y K A • R O B O T Y K A
3/2018 ISSN 1427-9126 Indeks 339512
Cena 25,00 zł w tym 5% VAT
Technical Sciences Quarterly | - .
!
5 15
Cena 25,00 zł w tym 5% VAT
Technical Sciences Quarterly | - .
!
W numerze:
3
4/2018 ISSN 1427-9126 Indeks 339512
W numerze:
Od Redakcji
3
Od Redakcji
> * # M ( '
5
) * ' ' # Q
0*€ +`
@ a
SposĂłb wyboru optymalnej trajektorii lotu manewru antykolizyjnego realizowanego w otoczeniu ruchomych przeszkĂłd
" < #
" ? C ) #
„ (
# '
'( #
23
# #$ () ' @ = @
33
Patryk Szywalski. Dawid Wajnert
41
; # @ $
1+ & (
' # * '
8 ( # ' < '
1 '
. U * ' ' ( U U 5&
1 # *
23
" $
25
?
31
4 # Â ' %
' 0U
0U # *
/ U *€ ' ' ' # '
'
( ~5. # ' # < =)# ? #
' * ' * ' '
1' ( ( v,
Ponadto:
+
Informacje dla Autorów – 69 | 1/// v! 0 3 <,+O3 2 #
) 4 ) 1/// $ 73 | E,
*$ 76 | $ * $ 80 | Nasze wydawnictwa – 84 | $ $ ! " " $
85 | Kalendarium – 86 | &," ( ) *$ 87 |
93 |
v 5 +,'- =% +,'- B €" 0 ) $ 0 0 +5 E5 94 |
1  +,'- = +,'-€ q EV," 0 ) $ 95 | 96
Informacje dla Autorów – 61 | PP-RAI 2018 – Pierwszy Zjazd Polskiego Porozumienia na Rzecz Rozwoju Sztucznej Inteligencji – 65 | Awans naukowy – ! " " $ # B ! 69 | Nasze wydawnictwa – 70 | Kalendarium – 71 | 17th 1q v B ! 0 v 00 20 72 | Konferencja Roboty Medyczne 2018 – 75 | Universities of the Future – 76 | # 1 $ ) +,'} 78
www.par.pl
www.automatykaonline.pl/automatyka
72
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
D; %D (K’“ K
0 % < % ! " ' ' " &' 3 $ 0 B $ $ $ B 3 $ # F $ )( )$ ) ! > !
# # ! 0# *$ $ ) 0 $ " ) 0 " ! " " ) q "
Prezentowana publikacja stanowi kompendium wiedzy z zakresu mechaniki punktu ) = # #
# M ) = Â? ) ^ M ]= Â? ! ^ ) ! ^ Â? # Â? ) $# # Â?
)
M ! ) mosfery ziemskiej i termomechaniki budowlanej. < ^H ]!
) ^ $ M &
< ) ] )& ?
Newtona przytoczono trzy prawa ruchu,
K % ]
$ ! ] ) Z ! ) ] # ][& (] +
$ ^
# #
&
] $ $# & ! )M ] $ =) ]! & # ) Z ] )
) # # ) #
$ [ A ) ) ] ) M # M A Oxyz oraz ,
# ] M ! &
] )& & = ' ) # M
^ = #
) = #$
= ' ) # # $# ] ) ]&
@ # ] ) #
# M ) A
^
) ] # #
# M ) ) =] $ H
#^ ^
^ ]= & ? ! ]
= ] # ) ] M
] &
M ! ! ] M !
modele. > A ) ) ) )
' ) $ )
# M ) M ) = ^ # # # = ] ) & )M # ) )
M! $ & * M $ ) ^
! ^ M
) ) ) ) $ -
^ & ' M ]
) # ) ) # # ] ) )& ? $ = )
$ ' ) A ) ) ] $ ) M!
= ' ) M !
=M
)
M ) & ) A
$ -
) # $ M$ = M ^ ] )
^ M ) & C#
)M
# # ) ) # )
M ) ) M '
Z D([ ) $ ^ '
Z ([& # ) ^
] ) ) )
M $ =
$ ^ M
M! ^ = ) ) M ] #!]
pomoc w praktycznych zastosowaniach M! M )
M & % ' # ! = # ) czono zestawienie pozycji literaturowych A #) UQG& =
) ^ = M
) # M M
M # M! M
! # & =]
H ! !
$ )# ] ) )
i obliczeniami komputerowymi z zakresu & )
M &
B ( " T Andrzej Flaga, & (
( '( , Wydawnictwo Naukowe PWN, 2019, ISBN 978-83-01-20179-1, str. 596.
# ] ] # M ' M! A
M ! )
$ # $ $ ) A ]
$ ' & retyczne analizy wymiarowej zilustrowano M
A )# =
) # $# ! !
#& ? $ ) ]=
^ M ) ] )
)
M A
! $
73
?˜* "<D%K j % W"(
B % 4 )!*+ młodzi młodzi
innowacyjni innowacyjni yj yj
= !
$ # & & ( ' * !& #
#
_D ) = $ =
`
) ! ^
= M D @
$ #
) = # ] ) # $
$ = & ?
] $ $ = # ! ( > M M M ]& # # )& & ! # ] # ] ) = E
]
] ]
H $ ^H #
) ^
M # ) ! ! H ) = $ ) # ] ) # ] ) ) = E
] ! ) = E
) =] $ H
) $ # M &
)
$ & * !& %
_
i# ž# X ž# X
) # = = )
)# # )
=
$ )`
#
$ ) )M )
( Z =& ( )# #
=[ # )
= M ] $# & #
$ ]
$ )# (
= # = ' M = # ]
= #$
! ^ # ] $ Z ) [
= ] ] & ?
# )# $ # ]
$# # M # ) #) !
= M ? i
= ) = ) ª"
M #) ^ # )
# ^ ) & = ) )
= ) ! ^ $ # ] '
) $ = $ &
?
) $ -
$(< " 0 6 6 ŒB % 4 ˜& ( % 0 4 8 0 < 0480 % % ( % R % S S
( % S F %
Ocenie konkursowej, jak co roku, poddano ^H
) Z) = ! [
+ # )
) M
$ =M ) # = )
M #
)
! '&
? ) # # # E hQ
A F
UU
) = OU
! & ?
M
@ ' A
H
) @M b# A )
& ) # $ #
$ M # > = #
Z > [ @ ZW [ @ Z \ ] [
Z \ [ Z [ #$ Z
#$ [ J Z JM [ ZW ? ) ' # [
Z [ Z ' [ ( Z< ) W C = [ ? Z ? ?
) C [ < @M
ZW < =M [& ; ) # ) ) )
# #+ • & $& !& # A K #
> ' ( ) % A
] • & !& ( A
& $
A ) K # # ) ) M K
TP
P
O
M
I
A
R
• !& = ' A
# =
) # ) " $ • & $& !& <$ = % A
K # > ' ( ) % • & & $& !& #
A ? ? • & & $& !& C A
? ( ) M " $ ? • & & $& !& ; < ' A
) # ) " $ % #
# )
= ! M!
Z
#) # = _ ) # ) " $ ` ) $ ! = _ WC C˜ ` # = _ #
# ) $ " $ © K = ( ) `[& ? $ ^ M = # #
OT ) OGUc &
ffKKK % # C _ # ) A % ^ ` WC CK % OGUc& ( )
$ M = ) ) ) ) $# # & > # #)
) ) ^ ^ ' ! M$
] ) ^
l & C )
$
# # $ )
# ) $ ) $ & < )
) )M
) ) = ] =
) ] ) ] =
Y
•
A
U
T
O
M
A
T
Y
K
A
•
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
) # ) " $
K((% UPOT cUOF "& Oh % UBOGUc Uczelnia
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
AGH Akademia GĂłrniczo-Hutnicza
2
5
5
3
19
4
4
3
5
3
7
Akademia Marynarki Wojennej
–
–
–
–
–
–
–
–
1
–
–
Akademia Morska w Gdyni Uniwersytet Morski w Gdyni
–
1
–
–
–
–
–
–
–
–
1
Akademia Techniczno-Humanistyczna w Bielsku-Białej
1
1
–
–
–
–
4
–
–
–
–
Instytut Badań Systemowych PAN
–
–
–
–
2
–
–
–
–
–
–
Państwowa Wyşsza Szkoła Zawodowa w Sanoku
1
–
–
–
–
–
–
–
–
–
–
Państwowa Wyşsza Szkoła Zawodowa w Tarnowie
–
–
–
1
2
–
2
–
–
–
–
Państwowa Wyşsza Szkoła Zawodowa w Zamościu
–
–
–
–
–
1
–
–
–
–
–
Politechnika Białostocka
–
1
–
3
–
1
–
–
5
6
2
Politechnika Częstochowska
–
–
–
–
1
–
–
–
–
–
–
Politechnika Gdańska
5
3
8
7
12
5
4
6
3
7
9
Politechnika Koszalińska
–
–
–
–
1
1
–
–
–
–
–
Politechnika Krakowska
–
5
6
1
3
–
–
1
–
–
1
Politechnika Lubelska
–
–
–
–
–
–
–
1
–
–
1
Politechnika Ĺ Ăłdzka
1
7
–
4
8
1
5
1
2
–
1
Politechnika Opolska
4
2
–
1
2
1
6
–
2
2
1
Politechnika Poznańska
2
2
3
–
10
1
6
2
5
1
1
Politechnika Rzeszowska
1
–
–
–
2
–
2
–
4
–
–
Politechnika ĹšlÄ…ska
7
5
6
4
10
8
2
6
5
7
1
Politechnika Świętokrzyska
–
–
–
1
2
–
–
–
–
–
1
Politechnika Warszawska
18
17
18
13
22
11
17
14
17
5
3
Politechnika Wrocławska
4
5
5
2
8
6
8
1
–
7
–
Polsko-Japońska Wyşsza Szkoła Technik Komputerowych
–
–
–
1
–
–
–
–
–
–
–
Uniwersytet Mikołaja Kopernika w Toruniu
–
1
–
–
1
1
–
1
1
–
–
Uniwersytet Rzeszowski
–
–
–
–
3
1
1
–
–
–
–
Uniwersytet Technologiczno-Przyrodniczy w Bydgoszczy
–
2
–
–
1
–
–
–
–
–
–
Uniwersytet Warmińsko-Mazurski w Olsztynie
–
1
–
–
–
–
–
–
–
–
1
Uniwersytet Warszawski
–
–
–
–
–
–
–
1
–
–
–
Uniwersytet ZielonogĂłrski
–
–
1
–
2
1
–
–
1
–
1
Wojskowa Akademia Techniczna
–
1
–
1
2
1
–
–
–
–
4
Wyşsza Szkoła Informatyki Stosowanej i Zarządzania w Warszawie
4
–
–
–
–
–
–
–
–
–
–
Zachodniopomorski Uniwersytet Technologiczny w Szczecinie
1
2
–
–
3
1
2
2
5
3
3
TÂ
%" hBOGUÂ
) M M
] ^H
= )M ) & = !& % $ # #
_ = $ # ) # ) = $ ) $ = `
$# ) )
= ) ] # ) = $ ) $ = &
# $
) $ = #) !
# # ) & ? )
$
= ] ) ) O* = $ #
# ) = $ ) $ = &
)
# ) ^ $ = $gKD?& @ M )
) # ) = $ ) $ = $ ) = $ #&
? )
$ M!
) & ? ! )
M! $ & ^)
= & >
M $
) $ ) ) H& ?
) # ] !
) ! $ # ]
) # ]
) ) ) ) & K ! @
# ] _" ] $
$ )M )
M
#
=
E ) ( C f`& ) $
=# ) # ' ) ] # #
# M
# ] ' # ) = E ) &
? # $M ) M
$ Z ) KC
C ) P&G ) = K # " [ ! ] '
) # #
# #=& ? $M )M ) '
) ) )&
Prace doktorskie
I Nagroda
dr inş. Michał Mikulski – Electromyography-based diagnostics and therapy of clinically weak patients Politechnika Śląska, Wydział Automatyki, Elektroniki i Informatyki Promotor – prof. dr hab. inş. Jarosław Śmieja
II Nagroda
dr inş. Michał Nowicki – Methods for the Fusion of quantitative and qualitative information using factor graph optimization for the simultaneous localization and mapping problem Politechnika Poznańska, Wydział Elektryczny Promotor – prof. dr hab. inş. Piotr Skrzypczyński
Wyróşnienie
dr inş. Anna Miller – Synteza elektronawigacyjnego układu sterowania automatycznego ruchem równoległym statków z wykorzystaniem metod predykcyjnych Uniwersytet Morski w Gdyni, Wydział Elektryczny Promotor – prof. dr hab. inş. Witold Gierusz
Wyróşnienie
dr inş. Paweł Latosiński – Ślizgowe sterowanie dyskretnymi obiektami dynamicznymi z wykorzystaniem predefiniowanych przebiegów zmiennej przełączających Politechnika Šódzka, Wydział Elektrotechniki, Elektroniki, Informatyki i Automatyki Promotor – prof. dr hab. inş. Andrzej Bartoszewicz Prace dyplomowe magisterskie
I Nagroda
mgr inş. Norbert Kukurowski – Mapowanie nieznanego obszaru z wykorzystaniem autonomicznego robota mobilnego Uniwersytet Zielonogórski, Wydział Informatyki, Elektrotechniki i Automatyki Promotor – prof. dr hab. inş. Krzysztof Patan
II Nagroda
mgr inş. Piotr Kwiatkowski – Zrobotyzowana mobilna stacja paletyzująca oparta na robocie współpracującym Politechnika Warszawska, Wydział Mechatroniki Promotor – prof. dr hab. inş. Mariusz Olszewski
Wyróşnienie
mgr inş. Szymon Grocholski – Rehabilitation Platform for Children with Cerebral Palsy Politechnika Gdańska, Wydział Elektroniki, Telekomunikacji i Informatyki Promotor – prof. dr hab. inş. Zdzisław Kowalczuk
Wyróşnienie
mgr inş. Andrzej Zakręcki – Opracowanie metody regeneracji łopatek turbiny parowej za pomocą technologii Laser Metal Deposition AGH Akademia Górniczo-Hutnicza, Wydział Inşynierii Mechanicznej i Robotyki Promotor – prof. dr hab. inş. Tadeusz Uhl
Wyróşnienie
mgr inş. Anna Legierska – Super-resolution algorithms of diffusion magnetic resonance imaging data AGH Akademia Górniczo-Hutnicza, Wydział Elektrotechniki, Automatyki, Informatyki i Inşynierii Biomedycznej Promotor – dr inş. Tomasz Pięciak Prace dyplomowe inşynierskie
I Nagroda
inş. Maciej Pilch – Badanie wpływu ciśnienia na charakterystykę emisji molekularnych sensorów luminescencyjnych w polimerowych materiałach powłokowych Politechnika Krakowska, Wydział Inşynierii i Technologii Chemicznej Promotor – dr inş. Joanna Ortyl, prof. dr hab. inş. Roman Popielarz
II Nagroda
inş. Paweł Gajewski – Rozwiązanie wybranych problemów informatycznych w projektach racjonalizacji technicznej procesów produkcyjnych przygotowania powierzchni prowadzonych w firmie SciTeex Politechnika Warszawska, Wydział Mechatroniki Promotor – prof. dr hab. inş. Mariusz Olszewski
Wyróşnienie
inş. Dawid Adamski – Wykonanie modelu laboratoryjnego układu sterowania i zobrazowania połoşenia satelity typu Cube Wojskowa Akademia Techniczna, Wydział Mechatroniki i Lotnictwa Promotor – mjr dr inş. Maciej Henzel
Wyróşnienie
inş. Tomasz Król – Projekt i realizacja rozproszonego systemu sterującego oświetleniem z analizą zuşycia energii i zdalnym zarządzaniem Politechnika Opolska, Wydział Elektrotechniki, Automatyki i Informatyki Promotor – dr inş. Grzegorz Bialic
Wyróşnienie
inş. Anna Wyşlińska – Prototyp system do monitorowania urazów głowy rowerzysty AGH Akademia Górniczo-Hutnicza, Wydział Elektrotechniki, Automatyki, Informatyki i Inşynierii Biomedycznej Promotor – dr inş. Andrzej Skalski
Wyróşnienie
inş. Maciej Kurcius – Opracowanie systemu sterowania dla 3-osiowej frezarki CNC AGH Akademia Górniczo-Hutnicza, Wydział Elektrotechniki, Automatyki, Informatyki i Inşynierii Biomedycznej Promotor – dr inş. Krzysztof Kołek
TF
P
O
M
I
A
R
Y
•
A
U
T
O
M
A
T
Y
K
A
•
% @ B ( " T
R
O
B
O
T
Y
K
A
NR 1 / 20 1 9
Nowe Studia Podyplomowe na Wydziale Mechatroniki
TRANSFORMACJA PRZEMYSŁOWA 4.0 Początek zajęć w październiku, w roku akademickim 2019/2020, zgłoszenia:
do 30 września 2019 Więcej informacji na stronie Rekrutacja na Studia na stronie
http://iair.mchtr.pw.edu.pl/studia https://rekrutacja.pw.edu.pl/pdp
CO GWARANTUJĄ? Wiedzę oraz kompetencje do rozwiązywania problemów przejścia od etapu zautomatyzowanej i zrobotyzowanej produkcji, charakteryzującej się liniową strukturą wytwarzania produktów, do etapu cyfryzacji i mechatronizacji produktów i produkcji o rozproszonej strukturze, produkcji zinformatyzowanej i zinternetyzowanej, na współczesnym i przewidywanym poziomie rozwoju tych technik, czyli przejścia do etapu Przemysłu 4.0. DLA KOGO? Dla absolwentów studiów wyższych I lub II stopnia na kierunkach technicznych, ekonomicznych lub menedżerskich. Wskazane, ale nie absolutnie obowiązujące, jest doświadczenie zawodowe z zakresu racjonalizacji technicznej lub zarządzania przemysłowymi procesami produkcyjnymi. Słuchaczem może być także osoba mająca profesjonalny kontakt z tymi dziedzinami, choć niekoniecznie formalnie wykształcona we wskazanych kierunkach studiów. DO KOGO KIEROWANE? Przede wszystkim do menedżerów zarówno małych, jak i dużych firm i ich oddziałów, bez względu na specyfikę produktową oraz do kierowników i pracowników zespołów projektujących cały cykl życia scyfryzowanych produktów lub wybranych etapów ich życia w łańcuchu PLM (Product Lifecycle Management), rozpoczynając od podania koncepcji i dokumentacji wirtualnej, decyzji o podjęciu produkcji, przejścia do środowiska realnego, opracowania wspomaganej komputerowo produkcji, montażu i logistyki, kontroli eksploatacji aż do sterowanego recyklingu.
KTO WSPOMAGA? Studia są wspomagane przez firmy przemysłowe o znaczeniu wykraczającym poza rynek europejski – są to holdingi Festo, Grupa Bosch-Rexroth, Fanuc, Balluf i Siemens. Wyraża się to w korzystaniu z ich zasobów programowych i sprzętowych oraz z udziału doświadczonych praktyków firmowych w prowadzonych na Studiach zajęciach wykładowych i laboratoryjnych. Planowane są także wyjazdy studialne, krajowe i zagraniczne, w tym 6-dniowy wyjazd studialny (Polska, Niemcy) oraz wyjazd na Leitmesse fuer Intelligente Automation und Robotik w Monachium. JAK PRZEBIEGAJĄ ZAJĘCIA? Jeden rok akademicki, dwa semestry, dziewięć zjazdów piątkowo-sobotnio-niedzielnych raz w miesiącu i dwa wyjazdy studialne – łącznie 250 godzin zajęć oraz praca końcowa przewidziana na 40 godzin. Zaliczenie studiów – na podstawie egzaminów z przedmiotów wykładowych i zaliczeń z pozostałych przedmiotów oraz obrony pracy końcowej. Dokumentem potwierdzającym ukończenie Studiów jest Świadectwo wydane przez Politechnikę Warszawską.
Adres i zgłoszenia: Instytut Automatyki i Robotyki na Wydziale Mechatroniki ul. św. A. Boboli 8, 02-525 Warszawa Elżbieta Boniecka, p. 253 tel. 22 234 85 55 fax 22 849 03 98 iair.podyplom@pw.edu.pl
GŁÓWNE TEMATY ZAJĘĆ: • Ewolucja technik wytwarzania produktów i środków ich produkcji • Platforma Przemysłu 4.0 • Automatyzacja i robotyzacja produkcji • Mechatronizacja, informatyzacja i internetyzacja produktu • Zadania zmechatronizowanego produktu w procesie produkcji i eksploatacji – Internet Rzeczy • Cyberbezpieczeństwo przemysłowych systemów informacyjnych • Standardy komunikacyjne automatyki przemysłowej • Polski przemysł przyszłości • Przykłady transformacji przemysłowej 4.0
43
* $ " ) ( @#
47
< =)# ? #
- ' ' ' 0 U
„
' * '
' 2 : 0*€ Y2 &
*