Карачун В.В., Мельник В.М., Саверченко В.Г.

Національний технічний університет України «КПІ»

ВИЗНАЧЕННЯ ПЕРЕДАТОЧНОЇ ФУНКЦІЇ ОБ'ЄКТА ПО РЕАКЦІЇ НА ЗАДАНИЙ ВХІДНИЙ СИГНАЛ

 

Для оцінки динамічних характеристик об’єкта, чи будь-якого іншого елемента автоматичної системи (який надалі будемо називати просто об’єктом), зручно використовувати фіксовану реакцію y(t) даного об’єкта на відомий (наперед заданий чи зафіксований синхронно з реакцією) вхідний сигнал x(t) .

Наближене значення передаточної функції вивчаємого каналу, як правило, шукають у вигляді.

                                      (1)

Їй відповідає диференціальне рівняння

                   (2)

Існує значна кількість способів визначення порядку n та коефіцієнтів передаточної функції (1), в тому числі і шляхом обробки сигналів x(t) та y(t), наприклад, метод площ Симою, метод послідовного інтегрування рівняння (2) за часом та інші.

Характерною особливістю розв’язуваної задачі є можлива присутність похибок у значеннях x(t) та y(t). Щоб нівелювати вплив цих похибок на результати ідентифікації, функції x(t) та y(t) згладжують, а це призводить до виникнення додаткових похибок згладжування, тобто коло проблем замикається.

Другою особливістю задачі є можлива невідповідність структури (1) глибинним особливостям динамічної поведінки досліджуваного об’єкту. Так, наприклад, порядок n нижче «справжнього». Остання обставина завжди має місце при апроксимації передаточною функцією типу (1) передаточних функцій більш високого порядку (задача пониження порядку систем). А для об’єктів з розподіленими параметрами «справжній» порядок взагалі дорівнює нескінченності.

Спроби ж «силового» приведення подібних об’єктів до структури (1) можуть мати наслідком проблеми з визначенням коефіціентів цієї структури (починаючи від поганої обумовленості відповідної системи алгебраічних рівнянь, яку доводиться розв’язувати, і закінчуючи «неприйнятними» значеннями коефіцієнтів, наприклад, такими, при яких об’єкт, що є явно стійким, може одержати знаменник своєї передаточної функції з коренями в правій півплощині комплексної площини, тощо).

При апроксимації кривих аналітичними залежностями, зокрема, при їх згладжуванні, широко використовується метод найменших квадратів (МНК), який несе в собі елемент згладжування. Отже, має сенс скористатися цією особливістю МНК для компенсації впливу похибок (в тому числі за рахунок невідповідності структур) при розв’язанні задачі ідентифікації.

Будемо вважати, що функції x(y) та y(t) задані масивами Х,У:CoefR їх ординат такої структури 

-1

    0

1

2

3

...

m

 

 

X

m

x0

    x1

          x2

x3

...

xm

 

 

 

-1

0

        1

2

3

...

m

 

 

Y

h

y0

y1

y2

y3

...

ym

 

 

 

Тут type CoefR = array[-1.. Mmax] of real;

      h - крок за часом, h = D/m;

      D - час спостереження сигналів x(t) та y(t);

      m - кількість кроків за часом.

За показник якості ідентифікації приймемо функцію

              (3)

Тут індекс і означає, що відповідна функція визначається в точці  t = ih.

Прирівнюючи до нуля похідні від Е0 відносно коефіцієнтів a1, a2,..., an, b0,b1,..., bn, одержимо систему рівнянь, розв’язком якої і будуть шукані коефіцієнти структури (1). Але ж, як ми вiдзначили вище, сигнали x(t) та y(t) можуть мати похибки, в тому числі високочастотні. А це може різко погіршувати точність визначення похідних, присутніх в формулі (3), особливо похідних високих порядків. Результатом буде погіршення якості ідентифікації.

Вплив високочастотних «шумів» можна спробувати ослабити шляхом k послідовних інтегрувань (0<k<n) рівняння (2) за часом t в межах від 0 до t.

Тоді рівняння (2) предстануть у вигляді

                           (4)

де

Тут використано позначення

За показник якості ідентифікації візьмемо

.                    (5)

Умови мінімізації Е

                                               (6)

Підставляючи значення Е з (5) в (6) та ділячи на 2 одержимо

           (7)

Якщо розташувати шукані коефіцієнти у такій послідовності a1, a2, ..., an, b0, b 1, b 2, ..., b n, то елементи розширеної матриці С системи (7) лінійних алгебраїчних рівнянь можна обчислювати за такими формулами

                   (8)

         (9)

Після того як масив С:Matr сформований (type Matr=array[1..2*Nmax+1, 1..2*Nmax+2] of real), викликаємо процедуру SystUr таким чином SystUr(2*n+1,C,A).

Масив А повертається у вигляді

 

-1

0

1

2

..

   n

n+1

n+2

...

2n+1

 

A

2n+1

~

1

2

..

an

  b0

  b1

...

  bn

 

 

 

Цей масив розбивається на два

 

-1

0

 1

 2

...

 n

 

 A

 n

1

a1

a2

...

an

 

 

 

-1

 0

1

 2

...

n

 

 B

 n

b0

b1

b2

...

bn

 

 

При потребі масиви А та В (особливо В) можна відредагувати, тобто скоригувати (в сторону зменшення) вміст комірки з номером  -1, коли старші коефіцієнти за модулем виявляться меншими деякої наперед заданої малої величини ε.

Наведемо описаний алгоритм у вигляді процедури WpXY. При цьому будемо вважати, що крок у часі h достатньо малий, щоб операцію інтегрування в формулах для e(t), f(t), g(t) рівняння (4) можна було б реалізувати за формулою трапецій, а диференціювання в тих же формулах - через ліві, праві та центральні різниці.

Procedure  WpXY(X,Y: CoefR; n,k:integer; var B,A:Coef);

     var  m,s,z,i     :integer;

          h,Sum:real;  M1,M2:CoefR; C:Matr;

     procedure Integr (var Jf:CoefR);

          var s:integer;

              Mf:CoefR;

          begin

              Mf:=Jf;  Jf[0]:=0;

              for s:=1 to m do

              Jf[s]:=Jf[s-1]+(Mf[s-1]+Mf[s])/2*h

            end;

     procedure Dif (var Df:CoefR);

          var s:integer;

              Mf:CoefR;

          begin

              Mf:=Df; Df[0]:=(Mf[1]-Mf[0])/h;

              Df[m]:=(Mf[m]-Mf[m-1])/h;

              for s:=1 to m-1 do    Df[s]:=(Mf[s+1]-Mf[s-1])/(2*h)

            end;

     procedure EzFz (z:integer; Mf:CoefR;  var Mz:CoefR);

          var i,j,f,s :integer;

              Iy:CoefR;

              t,Ts,R:real;

          begin

              for i:=0 to m do

                 begin

                     t:=i*h;

                     if z<=k-1   then    if   i=0

                         then R:=0     else

                             begin

                                 Iy:=Mf;

                                 for j:=1 to k-z do Integr(Iy);

                                 R:=Iy[i];

                                 for s:=0 to z-1 do

                                   begin

                                       Iy:=Mf;

                                       for j:=1 to s do  Dif(Iy);     Ts:=1;

                                       for j:=1 to k-z+s do

                                           Ts:=Ts*t/j;    R:=R-Iy[0]*Ts;

                                     end

                               end

                     else

                         begin

                             Iy:=Mf;

                             for j:=1 to z-k do Dif(Iy);     R:=Iy[i];

                             for s:=0 to k-1 do

                                begin

                                    Iy:=Mf;

                                    for j:=1 to z-k-s do Dif (Iy);

                                    Ts:=1;

                                    for j:=1 to s do

                                        Ts:=Ts*t/j;

                                    R:=R-Jy[0]*Ts

                                  end

                           end;

                         Mz[i]:=R

                   end

            end;

     Рrocedure Gt(var Mg:CoefR);

          var  s:integer;

          begin

               Mg:=Y;

               for s:=1 to k do Integr(Mg)

            end;

          begin

                m:=round(X[-1]); h:=Y[-1];

                for z:=1 to 2*n+1 do

                  for s:=1 to 2*n+2 do

                    begin

                       Sum:=0;

                       if z<=n         then

                            begin

                               EzFz(Z,Y,M1);

                               if s<=n      then EzFz(S,Y,M2)

                                  else          if  s<=2*n+1

                                        then  EzFz(s-n-1,X,M2)

                                        else  Gt(M2)

                              end

                         else      begin

                                  EzFz(z-n-1,X,M1);

                                  if s<=n   then  EzFz(S,Y,M2)

                                     else       if s<=2*n+1

                                          then  EzFz(s-n-1,X,M2)

                                          else  Gt(M2)

                               end;

                       for i:=0 to m do

                            Sum:=Sum+M1[i]*M2[i];

                       c[z,s]:=Sum

                      end;

                      Systur(2*n+1,C,A);

                      A[-1]:=n; A[0]:=1; B[-1]:=n;

                      for s:=0 to n do    B[s]:=A[n+1-s]

            end;