100724

Идентификация динамического объекта 2-го порядка с использованием метода наименьших квадратов

Лабораторная работа

Математика и математический анализ

Изучение особенностей идентификации математических моделей линейных динамических систем с использованием метода наименьших квадратов.

Русский

2018-02-10

267 KB

2 чел.

Идентификация динамического объекта 2-го порядка

с использованием метода наименьших квадратов

Цель работы: Изучение особенностей идентификации математических моделей линейных динамических систем с использованием метода наименьших квадратов.

По варианту задана следующая передаточная функция:

1. Для динамического объекта, заданного передаточной функцией с параметрами, согласно варианту задания, используя пакетMatLabполучить выходной сигнал при входном воздействии в виде единичной ступенчатой функции (период дискретизацииΔt = 0,1 с ).

2. Получить реализацию выходного сигнала объекта при наличии шума измерения (амплитуда шума sv = 0.1).

Программа для реализации пунктов 1 и 2:

clear all

clc

K0 = 20;

T1 = 2.8;

T1kv = T1^2;

T2 = 13;

W1 = tf(K0, [T1kv T2 1]);

Tend = 60;

dt = 0.1;

t = 0 : dt : Tend;

N = length(t);

u = ones(N,1);

y0 = lsim(W1,u,t);

sv = 0.1;

v = sv*randn(N,1);

ys = y0 + v;

3. Используя функцию plot построить в одной системе координат графики переходного процесса без учёта и с учётом шума измерений. Результат представлен на рис.1.

Программа для реализации пункта 3:

figure (1);

hold on;

plot (t, y0);

plot (t,ys);

grid;

hold off;

Рисунок 1 - Графики переходных процессов системы без учёта и с учётом шума измерений

4. Разработать и сохранить в текущей директории m-файл функцию, реализующую регрессионный МНК для определения параметров дискретной модели заданного динамического объекта и соответствующих им параметров аналоговой модели.

Разработанная функция для регрессионного МНК:

function [ K0m, T1kvm, T2m ] = F_MNKreg( y,u,N,dt )

S1 =sum(y(1:N).^2);

S2 = sum(y(2:N).*y(1:N-1));

S3 = sum(y(2:N-1).^2);

S4 = sum(y(3:N).*y(1:N-2));

S5 = sum(y(2:N).*u(1:N-1));

S6 = sum(y(2:N-1).*y(1:N-2));

S7 = sum(y(1:N-1).*u(1:N-1));

S8 = sum(y(1:N-2).^2);

S9 = sum(y(1:N-2).*u(2:N-1));

S10 = sum(u(1:N-1).^2);

A = [S3, S6, S7; S6, S8, S9; S7, S9, S10];

B = [S2; S4; S5];

C = A^(-1) * B;

T1kvm = (dt^2) / (1 - C(1) - C(2));

T2m = ((C(2)+1) * T1kvm + dt^2) / dt;

K0m = C(3) * T1kvm / dt^2;

end

Программа для расчета параметров модели объекта:

[K0m0,T1kvm0,T2m0] = F_MNKreg(y0,u,N,dt);

[K0ms,T1kvms,T2ms] = F_MNKreg(ys,u,N,dt);

disp ('Параметры с интенсивностью шума 0.1');

disp ('K0m0,         T1kvm0,      T2m0');

disp ([K0m0,T1kvm0,T2m0]);

disp ('K0ms,         T1kvms,      T2ms');

disp ([K0ms,T1kvms,T2ms]);

Использование данной программы вMatlab:

K0m0,         T1kvm0,      T2m0

   19.9747   11.0015   14.1726

K0ms,         T1kvms,      T2ms

   20.2951    1.0156   15.1819

5.Определить амплитуду шума, при которой оценка параметров ухудшается примерно на 15%, и в дальнейшем использовать данное значение амплитуды шума при выполнении вычислительных экспериментов.

Получим реализации выходного сигнала и построим в имеющейся системе координат графики переходного процесса для полученных оценок параметров объекта без учёта и с учётом шума измерений. Используя обращение к разработанной m-файл функции получим оценку параметров объекта идентификации с учётом и без учёта шума измерений(рис.2), с помощью следующей программы:

sv = 0.5;

v = sv*randn(N,1);

ys = y0 + v;

figure (2);

hold on;

plot (t, y0);

plot (t,ys,'r');

grid;

hold off;

Рисунок 2 - Графики переходных процессов системы без учёта и с учётом шума измерений

6. Рассчитать  и вывести на экран математическое ожидание и среднеквадратичное значение ошибки выходного сигнала объекта с параметрами, полученными в результате идентификации с учетом и без учета шума измерений. Программа:

W2 =tf(K0m0,[T1kvm0T2m0 1]);

y2 =lsim(W2,u,t);

M2 = mean(y2-y0);

sx2 = std(y2-y0);

disp('Математическое ожидание');

disp(M2);

disp('Среднеквадратичное отклонение');

disp(sx2);

W2s = tf(K0ms,[T1kvms T2ms 1]);

y2s = lsim(W2s,u,t);

M2s = mean(y2s-y0);

sx2s = std(y2s-y0);

disp('Математическое ожидание');

disp(M2s);

disp('Среднеквадратичное отклонение');

disp(sx2s);

Результат выполнения программы:

Математическое ожидание

   -0.0446

Среднеквадратичное отклонение

    0.0427

Математическое ожидание

    0.2440

Среднеквадратичное отклонение

    0.7983

7. Произвести сглаживание зашумлённой переходной характеристики согласно методике, рассмотренной в п.5 теоретической части.

m = 13;

y = zeros(size(t));

y(1) = ys(1);

y(N) = ys(N);

for j = 2:m

    y(j) = sum(ys((1):(j*2-1)))./(2*j-1);

end

for j = m+1:N-m

    y(j) = sum(ys((j-m):(j+m)))./(2*m+1);

end

for j = N-m+1:N-1

    y(j) = sum(ys((j*2-N):(N)))./(2*(N-j)+1);

end

y = y';

figure (3);

plot (t,y)

Рисунок 3 - График сглаживания зашумленной переходной характеристики

8. Используя обращение к разработанной ранее m-файл функции получить оценку параметров объекта идентификации после сглаживания шума измерений. Вывести на экран полученные параметры передаточной функции. Получить реализацию выходного сигнала и построить в имеющейся системе координат график переходного процесса. Рассчитать и вывести на экран математическое ожидание и среднеквадратическое значение ошибки выходного сигнала объекта с параметрами, полученными в результате идентификации после сглаживания шума измерений.

Программа  и график переходного процесса:

[K0m01,T1kvm01,T2m01] = F_MNKreg(y,u,N,dt);

disp('Оценка параметров после сглаживания');

disp ('K, T1, T2');

disp ([K0m01,T1kvm01,T2m01]);

W3 = tf(K0m01,[T1kvm01 T2m01 1]);

y3= lsim(W3,u,t);

M3 = mean(y3-y0);

sx3 =std(y3-y0);

disp('Математическое ожидание');

disp(M3);

disp('Среднеквадратичное отклонение');

disp(sx3);

figure (4);

plot (t, y3);

grid;

Результат выполнения программы:

Оценка параметров после сглаживания

K, T1, T2

   20.0602    2.0928   14.3697

Математическое ожидание

   -0.0274

Среднеквадратичное отклонение

    0.1618

Рисунок 4 - График сглаженной переходной характеристики

9. Разработать и сохранить в текущей директории m-файл функцию, реализующую явный МНК, описанный в разделе 4 теоретической части.

Программа:

function [K0m,T1kvm,T2m] =F_MNKnat(y,u,N,dt)

C = y(3:N);

B = zeros(N-2,3);

for i = 1:N-2

    B(i,1) = y(i+1);

end

for i = 1:N-2

    B(i,2) = y(i);

end

for i = 1:N-2

    B(i,3) = u(i+1);

end

A = (inv(B'*B))*B'*C;

a1 = A(1);

a2 = A(2);

b = A(3);

T1kvm = (dt^2)/(1-a1-a2);

T2m = ((a2+1)*T1kvm+dt^2)/dt;

K0m = (b*T1kvm)/(dt^2);

end

10. Для оценок параметров передаточной функции, полученных на основе применения явного МНК выполнить пп. 5 – 8 задания. Переходные процессы изображены на рис.5-6.

sv = 0.5; %

v =sv*randn(N,1);

ys = y0 + v;

[K0ms_1,T1kvms_1,T2ms_1] =F_MNKreg(ys,u,N,dt)

disp ('Параметры с интенсивностью шума 0.5')

disp ('K0m0,         T1kvm0,      T2m0')

disp ([K0m0,T1kvm0,T2m0])

disp ('K0ms,         T1kvms,      T2ms')

disp ([K0ms,T1kvms,T2ms])

Wn_s1 = tf(K0ms_n,[T1kvms_n T2ms_n 1])

yn_s1 = lsim(Wn_s1,u,t)

Mn_0 = mean(yn_0-y0)

sxn_0 = std(yn_0-y0)

disp('Математическое ожидание')

disp(Mn_0)

disp('Среднеквадратичное отклонение')

disp(sxn_0)

Mn_s1 = mean(yn_s1-y0)

sxn_s1 = std(yn_s1-y0)

disp('Математическое ожидание')

disp(Mn_s1)

disp('Среднеквадратичное отклонение')

disp(sxn_s1)

Результат выполнения кода:

Математическое ожидание

   -0.0324

Среднеквадратичное отклонение

0.0228

Математическое ожидание

    0.2133

Среднеквадратичное отклонение

0.8417

m = 13;

y_nat = zeros(size(t));

y_nat(1) = yn_s1(1);

y_nat(N) = yn_s1(N);

for j = 2:m

    y_nat(j) = sum(yn_s1((1):(j*2-1)))./(2*j-1);

end

for j = m+1:N-m

    y_nat(j) = sum(yn_s1((j-m):(j+m)))./(2*m+1);

end

for j = N-m+1:N-1

    y_nat(j) = sum(yn_s1((j*2-N):(N)))./(2*(N-j)+1);

end

y_nat = y_nat';

figure (6);

plot (t,y_nat);

Рисунок 5 - График сглаживания зашумленной переходной характеристики

[K0m0_nat,T1kvm0_nat,T2m0_nat] = F_MNKnat(y_nat,u,N,dt);

disp('parameters after smoothing');

disp ('K0m0_1,         T1kvm0_1,      T2m0_1');

disp ([K0m0_nat,T1kvm0_nat,T2m0_nat]);

W_nat = tf(K0m0_nat,[T1kvm0_nat T2m0_nat 1]);

y0_nat = lsim(W_nat,u,t);

M_nat = mean(y0_nat-y0);

sx_nat = std(y0_nat-y0);

disp('mean =');

disp(M_nat);

disp('std =');

disp(sx_nat);

figure (7);

plot (t, y0_nat);

Рисунок 4 - График сглаженной переходной характеристики

Вывод: изучили особенности применения метода наименьших квадратов для идентификации математических моделей линейных динамических систем на примере передаточной функции второго порядка, представленной в исходных данных.

При выполнении поставленной задачи было выявлено, что при увеличении интенсивности шума портятся параметры модели. Операция сглаживания приближает изменившиеся параметры к исходным, что видно из результатов вычисления мат.ожидания и среднеквадратичного отклонения. Из результатов работы видно, что явный МНК точнее регрессионного, так как математическое ожидание и среднеквадратичное отклонение при использовании данного метода без сглаживания точнее.


 

А также другие работы, которые могут Вас заинтересовать

69892. Дослідження ефективності алгоритму 3.19 MB
  Мета: набуття навичок визначення часової складності алгоритму. Теоретичні питання план Функція складності алгоритму. Види функції складності алгоритмів. Часова функція складності.
69893. Алгоритми сортування даних в оперативній пам’яті 59.5 KB
  Під сортуванням розуміють процес перестановки об’єктів даної множини в певному порядку. Мета сортування – полегшити подальший пошук елементів у відсортованій множині. В цьому значенні сортування присутнє майже у всіх задачах обробки інформації.
69894. Програмування задач обробки структур даних, розташованих на зовнішніх носіях 94 KB
  Мета: ознайомитися з поняттям файлу навчитися створювати і читати файли. Теоретичні питання план Поняття файлу. Основні типи файлових структур. Особливості роботи з текстовими файлами.
69895. Циклические вычислительные процессы 145 KB
  Вычислить значения функции a=1.6x3-1.5 на интервале (-1,1) с шагом изменения аргумента 0.25. Выдать на печать отрицательные значения функции с соответствующими им значениями аргумента. Вычислить и вывести на печать таблицу значений функции
69896. Создание псевдонима базы данных 1.11 MB
  При работе с таблицами локальных БД или СУБД сама база размещается либо в каталоге на диске и хранится в виде отдельного набора файлов, либо на удаленном сервере. Обращение к базе данных происходит по ее псевдониму (Database Alias).
69897. ЭЛЕКТРОННАЯ ТАБЛИЦА EXCEL. ЭЛЕМЕНТЫ ОКНА.НАСТРОЙКА ПАНЕЛИ ИНСТРУМЕНТОВ. ВВОД И РЕДАКТИРОВАНИЕ ДАННЫХ 465.5 KB
  Создайте электронные таблицы для расчета оплаты труда рабочих строительной бригады пропорционально отработанному времени, разряду и тарифам. Данные для расчетов разместить на разных листах.
69898. Термометры расширения и термометры манометрические 281 KB
  Самые старые устройства для измерения температуры жидкостные стеклянные термометры используют термометрическое свойство теплового расширения тел. Основные достоинства стеклянных жидкостных термометров простота употребления и достаточно высокая точность измерения...
69899. Безусловная одномерная оптимизация 503.5 KB
  Цель работы: знакомство с оптимизационными задачами, изучение различных методов одномерной оптимизации и сравнение эффективности их применения для конкретных целевых функций.
69900. Переменные. Операторы. Массивы 26.37 KB
  Цель работы: Создание простого приложения на Java, выполняющего небольшие вычисления с выводом результатов на консоль. Лабораторная работа 1: Создайте класс с названием Calc и метод main() в нем. Создайте в методе main() локальную переменную i типа int.