Ідентифікація

 

Продовження таблиці 3.1

Номер точки	Сигнал  Y(t)		Номер точки	Сигнал  Y(t)
	t	Y(t)		t	Y(t)
91	130,225724	1,39381767	98	138,4322268	1,126083919
92	131,167454	0,889482	99	139,2394237	1,150989385
93	132,37825	1,22881896	100	140,4502191	0,789860135
94	133,589045	1,36268584	101	141,5264817	1,11674437
95	134,79984	0,95797203	102	143,1408756	1,306648544
96	136,010636	0,8708029	103	144,0826054	1,026462057
97	137,086899	1,4809868	104	145,158868	1,070046622

 

3.2 Побудова моделі за допомогою поліноміальної інтерполяції поліномом 7-го степеня

Побудова моделі за допомогою  поліноміальної інтерполяції поліномом 7-го степеня виконуємо тим же способом як і для вхідного сигналу U(t), що описаний в першому розділі.

Нижче наведено лістинг програмного  модуля в середовищі Matlab:

 

clear;

load wt3.mat

for i=1:8;

    for j=1:8;

        A(i,j)=sum(x.^(i+j-2));

    end;

end;

for i=1:8;

    B(i)=sum(x.^(i-1).*y);

end;

a=B/A;

t=1:150;

f=a(8)*t.^7+a(7)*t.^6+a(6)*t.^5+a(5)*t.^4+a(4)*t.^3+a(3)*t.^2+a(2)*t+a(1);

hPlot=plot(t,f);

hold on;

plot(x,y,'red')

xlabel('t')

ylabel('Y(t)')

axis([0,150,-0.1,2.8]);

set( hPlot, 'LineWidth', 2 );

 

Нижче наведемо графічну модель (рис. 3.1):

Рис. 3.1. Поліноміальна модель вихідного сигналу Y(t) поліномом n-го порядку

 

 

Розв’язавши систему, отримаємо  такі значення коефіцієнтів:

 

a₀ = 1,02538836673611;

a₁ = 0,0268564689044343;

a₂ = -0,00252492654257438;

a₃= 9,12939501281913e-05;

a₄ = -1,58612302537886e-06;

a₅ = 1,41286442406615e-08; 

a₆ = -6,21079558699196e-11;

a₇= 1,06636334655407e-13.

 

3.3 Побудова моделі за допомогою полінома Чебишева

Побудова моделі за допомогою  полінома Чебишева виконуємо тим же способом як і для вхідного сигналу U(t), що описаний в першому розділі.

Нижче наведено лістинг програмного  модуля в середовищі Matlab:

 

clear;

load wt3.mat;

Cheb_por=7;

n=104;

fi=ones(n,7);

for i=1:n

    fi(i,2)=x(i)-(1/n)*sum(x);

end;

for p=3:Cheb_por; i=1:n;

        beta(p)=-sum(x(i).*fi(i,p-1).^2)/sum(fi(i,p-1).^2);

        gama(p)=-sum(x(i).*fi(i,p-1).*fi(i,p-2))/sum(fi(i,p-1).^2);

        fi(i,p)=(x(i)+beta(p)).*fi(i,p-1)+gama(p).*fi(i,p-2);

end;

A=fi'*fi;

B=fi'*y;

a=A\B;

for i=1:n; j=1:7; 

        Cheb(i)=a(j)'*fi(i,j)';

end;

hPlot=plot(x,Cheb);

hold on;

plot(x,y,'red');

xlabel('t');

ylabel('Y(t)');

set( hPlot, 'LineWidth', 2 );

 

Нижче наведено модель досліджуваного об’єкта (рис. 3.2).

Розв’язавши систему, отримаємо  такі значення коефіцієнтів:

a₁ = 1,05675426115593;           a₅ = 4,21939951118780e-08;

a₂ = -0,00317011396636858;           a₆ = -3,81381034627028e-10;

a₃ = -6,21837031383335e-05;           a₇ = -5,66801060511657e-12.

a₄ = 2,40061957256621e-06;

 

Рис. 3.2. Модель вихідного сигналу Y(t), побудована поліномом Чебишева

 

3.4 Побудова моделі за допомогою перетворення Фур’є

Побудова моделі за допомогою  перетворень Фур’є виконуємо  тим же способом як і для вхідного сигналу U(t), що описаний в першому розділі.

Нижче наведено лістинг програмного  модуля в середовищі Matlab:

 

clear;

load wt3.mat;

n=104;

Fur_por=4;

m=Fur_por;

for k=1:m+1

    A=(1/n)*sum(y);

    B(k)=(2/n)*sum(y.*cos(k*x));

    C(k)=(2/n)*sum(y.*sin(k*x));

end

for k=1:m+1

    for i=1:n

        FS(i,1)=A+sum(B(k)*sin(k*x(i))+C(k)*cos(k*x(i)));

    end

end

plot(x,y,'red')

hold on;

hPlot=plot(x,FS);

xlabel('t')

ylabel('Y(t)')

axis([0,150,-0.1,2.1]);

set( hPlot, 'LineWidth', 2 );

 

Розв’язавши систему, отримаємо  такі значення коефіцієнтів:

 

a₀ = 1,06570612602025;

 

 

Побудована модель представлена нижче, на рис. 3.3.

Рис. 3.3. Модель вихідного сигналу Y(t), побудована за допомогою перетворення  Фур’є

 

3.5 Статистична обробка даних

Розрахуємо коефіцієнти  дисперсії, кореляції та коваріації.

Розрахунок виконуємо  тим же способом як і для вхідного сигналу U(t), що описаний в першому розділі.

Нижче наведено лістинг програмного  модуля в середовищі Matlab:

 

clear;

load wt3.mat;

n=104;

i=1:n;

tser=sum(x)/n;

fser=sum(y)/n;

cov=sum((x-tser).*(y-fser))/(n-1)

Dt=sum((x-tser).^2)/(n-1)

Du=sum((y-fser).^2)/(n-1)

r=cov/sqrt(Dt*Du)

 

Після виконаного розрахунку отримаємо:

 

cov = -0.0393;

Dt = 1.7398e+003;

Du = 0.0574;

r = -0.0039

3.6 Знаходження періодограми сигналів

Побудову періодограми виконуємо тим же способом як і для вхідного сигналу U(t), що описаний в першому розділі.

Нижче наведено лістинг програмного  модуля в середовищі Matlab:

 

clear;

load wt3.mat;

n=104;

for k=1:n

    w(k,1)=(2*pi*(k))/n;

    U(k,1)=(1/sqrt(n))*sum(y.*exp(-j*w(k).*x));

end

Ui=imag(U);

Ur=real(U);

subplot(1,2,1);

plot(w,Ui);

title('a)');

xlabel('Нормована частота (Гц)');

ylabel('Спектр щільності (Дб/Гц)');

grid on

subplot(1,2,2);

plot(w,Ur);

title('б)');

xlabel('Нрмована частота (Гц)');

ylabel('Спектр щільності (Дб/Гц)');

grid on

 

 

Рис. 3.4 - Періодограма вихідного сигналу Y(t), де а) – уявна частина,

б) – дійсна частина 

3.7 Вибір оптимальної моделі

Вибір оптимальної моделі виконуємо тим же способом як і  для вхідного сигналу U(t), що описаний в першому розділі.

 

Оптимальні значення моделей

Таблиця 3.2

Вид функції	Поліном	Поліном Чебишева	Функція Фур’є
Значення критерію	0.0085	4.2762e-015	0.0639

Нижче наведено лістинг програмного  модуля в середовищі Matlab:

 

for i=1:a;n=104;

    Spol_Yt=sqrt(sum((y(i)-f(i)).^2/(n-2)))

end;

     SCheb_Yt=sqrt((sum(y-Cheb').^2)/(n-2))

 

     SFur_Yt=sqrt((sum(y-FS).^2)/(n-2))

 

Як видно з таблиці  3.2, значення критерію мінімальне у тому випадку, коли апроксимуючою функцією виступає Поліном Чебишева, отже модель побудована даним чином буде оптимальною.

 

Розділ  4. Побудова загальної моделі системи

 

4.1 Табулювання двох вхідних та одного вихідного сигналів

Табулювання графіків функцій (сигналів) наведені в попередніх розділах, тому таблиці та функції табулювання  вказувати не будемо. Лише зазначимо  деякі особливості.

Розіб‘ємо часову вісь t на N проміжків з інтервалом 2 секунди. Для кожного інтервалу визначимо значення у_р оптимальних моделей для сигналів U(t), V(t), Y(t).

        i:=1.. n       t_i:=2^.i

 

для U(t) – оптимальною моделлю є модель побудована за допомогою полінома Чебишева;

V(t) – оптимальною моделлю є модель побудована за допомогою полінома 7-ї степені;

Y(t) – оптимальною моделлю є модель побудована за допомогою полінома Чебишева.

        Спробуємо  визначити, яка з слідуючих моделей буде оптимальною для даної системи:

 

1. Y = a₀+a₁ U(t)+a₂W(t);
2. Y = a₀+a₁ U(t)+a₂W(t)+ a₃U(t)×W(t);
3. Y = a₀ +b ^.cos(U(t)) + d ^.sin(W(t)).

Для цього використаємо метод  найменших квадратів.

 

4.2 Побудова першої моделі системи, виду

Y = a₀+a₁ U(t)+a₂W(t)

За методом найменших квадратів  система рівнянь для першої моделі буде мати вигляд:

 

 

Нижче приведемо лістинг модуля рішення даної системи:

 

clear;

load wt1.mat

for i=1:8;

    for j=1:8;

        A(i,j)=sum(x.^(i+j-2));

    end;

end;

for i=1:8;

    B(i)=sum(x.^(i-1).*y);

end;

a1=B/A;

t=1:150;

f1=a1(8)*t.^7+a1(7)*t.^6+a1(6)*t.^5+a1(5)*t.^4+a1(4)*t.^3+a1(3)*t.^2+a1(2)*t+a1(1);

load wt2.mat

for i=1:8;

    for j=1:8;

        A(i,j)=sum(x.^(i+j-2));

    end;

end;

for i=1:8;

    B(i)=sum(x.^(i-1).*y);

end;

a2=B/A;

f2=a2(8)*t.^7+a2(7)*t.^6+a2(6)*t.^5+a2(5)*t.^4+a2(4)*t.^3+a2(3)*t.^2+a2(2)*t+a2(1);

load wt3.mat

for i=1:8;

    for j=1:8;

        A(i,j)=sum(x.^(i+j-2));

    end;

end;

for i=1:8;

    B(i)=sum(x.^(i-1).*y);

end;

a3=B/A;

f3=a3(8)*t.^7+a3(7)*t.^6+a3(6)*t.^5+a3(5)*t.^4+a3(4)*t.^3+a3(3)*t.^2+a3(2)*t+a3(1);

sf1=0;sf2=0;sf3=0;sf12=0;sf22=0;sf1f2=0;sf1f3=0;sf2f3=0;

n=150;

sf1=sum(f1);

sf2=sum(f2);

sf3=sum(f3);

sf12=sum(f1.^2);

sf22=sum(f2.^2);

sf1f2=sum(f1.*f2);

sf1f3=sum(f1.*f3);

sf2f3=sum(f2.*f3);

A1=[n sf1 sf2;sf1 sf12 sf1f2;sf2 sf1f2 sf22];

B1=[sf3 sf1f3 sf2f3];

a=B1/A1;

sysmod1=a(3).*f2+a(2).*f1+a(1);

hPlot=plot(t,f3,'red');

hold on;

hPlot1=plot(t,sysmod1);

grid on;

xlabel('t')

ylabel('Y(t)')

set( hPlot,'LineWidth', 2 );

set( hPlot1,'LineWidth', 2 );

S1=sqrt((sum(f3-sysmod1)^2)/(n-2));

 

Нижче побудуємо її графічне представлення  і знайдемо коефіцієнти а₀, а₁ , і а₂ , які визначатимуть модель:

 

a₀ =0,0288517639895134;

a₁ =0,567158190502034;

a₂ =0,424033755053556.

 

Рис. 4.1 – Загальна модель системи виду

Y = a₀+a₁ U(t)+a₂W(t)

 

4.3 Побудова другої моделі системи, виду

Y = a₀+a₁ U(t)+a₂W(t)+a₃U(t)W(t)

За методом найменших квадратів  система рівнянь для першої моделі буде мати вигляд:

 

 

Нижче приведемо лістинг модуля рішення даної системи:

 

clear;

load wt1.mat

for i=1:8;

    for j=1:8;

        A(i,j)=sum(x.^(i+j-2));

    end;

end;

for i=1:8;

    B(i)=sum(x.^(i-1).*y);

end;

a1=B/A;

t=1:150;

f1=a1(8)*t.^7+a1(7)*t.^6+a1(6)*t.^5+a1(5)*t.^4+a1(4)*t.^3+a1(3)*t.^2+a1(2)*t+a1(1);

load wt2.mat

for i=1:8;

    for j=1:8;

        A(i,j)=sum(x.^(i+j-2));

    end;

end;

for i=1:8;

    B(i)=sum(x.^(i-1).*y);

end;

a2=B/A;

f2=a2(8)*t.^7+a2(7)*t.^6+a2(6)*t.^5+a2(5)*t.^4+a2(4)*t.^3+a2(3)*t.^2+a2(2)*t+a2(1);

load wt3.mat

for i=1:8;

    for j=1:8;

        A(i,j)=sum(x.^(i+j-2));

    end;

end;

for i=1:8;

    B(i)=sum(x.^(i-1).*y);

end;

a3=B/A;

f3=a3(8)*t.^7+a3(7)*t.^6+a3(6)*t.^5+a3(5)*t.^4+a3(4)*t.^3+a3(3)*t.^2+a3(2)*t+a3(1);

n=150;

hold off;

sf1=sum(log(f1));

sf2=sum(log(f2));

sf1f1=sum(log(f1).^2);

sf1f2=sum(log(f1).*log(f2));

sf1sf2=sum(sf1.*sf2);

sf1f1f2=sum((log(f1).^2).*log(f2));

sf2f2=sum(log(f2).^2);

sf2f2f1=sum((log(f2).^2).*log(f1));

sf2f2sf1f1=sum((log(f1).^2).*(log(f2).^2));

sf3=sum(log(f3));

sf1f3=sum(log(f1).*log(f3));

sf2f3=sum(log(f2).*log(f3));

sf3f3=sum(log(f3).^2);

A1=[n sf1 sf2 sf1sf2;sf1 sf1f1 sf1f2 sf1f3;sf2 sf1f2 sf2f2 sf2f3;sf1sf2 sf1f1f2 sf2f2f1 sf2f2sf1f1];

B1=[sf3 sf1f3 sf2f3 sf3f3];

a=B1/A1;

sysmod2=a(4).*f1.*f2+a(3).*f2+a(2).*f1+a(1);

hPlot=plot(t,f3,'red');

hold on;

hPlot1=plot(t,sysmod2);

grid on;

xlabel('t')

ylabel('Y(t)')

set( hPlot,'LineWidth', 2 );

set( hPlot1,'LineWidth', 2 );

S2=sqrt((sum(f3-sysmod2)^2)/(n-2));

 

Розв’язавши її, знайдемо коефіцієнти а₀, а₁, а₂ і а₃, які визначатимуть модель:

 

a₀ = 0,00317016777990134;

a₁ =0,665835193742437;

a₂ =0,583053398632476;

a₃ =0,0181216832221273.

 

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

 

 

Рис. 4.2 – Загальна модель системи виду

Y = a₀+a₁ U(t)+a₂W(t)+a₃W(t)W(t)

 

4.4 Побудова третьої моделі системи, виду

Y = a₀ +b×cos(U(t)) + d×sin(W(t))

За допомогою методу найменших  квадратів знаходимо систему  рівнянь, розв’язавши яку можна  отримати коефіцієнти а₀, а₁, та а₂. За допомогою знайдених коефіцієнтів можна побудувати третю модель.

Система рівнянь буде мати вигляд:

 

Нижче приведемо лістинг модуля рішення даної системи:

 

clear;

load wt1.mat

for i=1:8;

    for j=1:8;

        A(i,j)=sum(x.^(i+j-2));

    end;

end;

for i=1:8;

    B(i)=sum(x.^(i-1).*y);

end;

a1=B/A;

t=1:150;

f1=a1(8)*t.^7+a1(7)*t.^6+a1(6)*t.^5+a1(5)*t.^4+a1(4)*t.^3+a1(3)*t.^2+a1(2)*t+a1(1);

load wt2.mat

for i=1:8;

    for j=1:8;

        A(i,j)=sum(x.^(i+j-2));

    end;

end;

for i=1:8;

    B(i)=sum(x.^(i-1).*y);

end;

a2=B/A;

f2=a2(8)*t.^7+a2(7)*t.^6+a2(6)*t.^5+a2(5)*t.^4+a2(4)*t.^3+a2(3)*t.^2+a2(2)*t+a2(1);

load wt3.mat

for i=1:8;

    for j=1:8;

        A(i,j)=sum(x.^(i+j-2));

    end;

end;

for i=1:8;

    B(i)=sum(x.^(i-1).*y);

end;

a3=B/A;

f3=a3(8)*t.^7+a3(7)*t.^6+a3(6)*t.^5+a3(5)*t.^4+a3(4)*t.^3+a3(3)*t.^2+a3(2)*t+a3(1);

n=150;

Fur_por=3;

m=Fur_por;

for k=1:m+1

    Am=(1/n)*sum(f3);