Практическое задание: Идентификация динамических систем по методу МНК

Идентификация динамических систем по методу МНК

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

ТЮМЕНСКИЙ ГОСУДАРСТВЕННЫЙ НЕФТЕГАЗОВЫЙ УНИВЕРСИТЕТ

ИНСТИТУТ ГЕОЛОГИИ И НЕФТЕГАЗОДОБЫЧИ

Кафедра "Автоматизации и управления"

ОТЧЕТ

по дисциплине

"Идентификация и диагностика систем"

к лабораторной работе № 3

"ИДЕНТИФИКАЦИЯ ДИНАМИЧЕСКИХ СИСТЕМ ПО МЕТОДУ МНК"

Тюмень 2015

Задание

Задана система из двух линейных дифференциальных уравнений:

где ,  - переменные состояния,  и  - входные сигналы,  и  - известные параметры.

Значения параметров  и  и выбор переменной состояния доступной к измерению указаны в Таблице 3.1 Параметры моделирования (временной интервал и шаг дискретизации) и форма входных сигналов  и  указаны в Таблице 3.2.

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

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

Значения параметров системы

Таблица 3.1

Параметры моделирования и формы сигналов управления

Таблица 3.2

Выполнение

Составим схему численного моделирования системы

Рассмотрим модель в форме пространства состояний:

Здесь  - переменные состояния,  - сигналы управления,  и  - параметры модели.

Зададим временной период  и шаг интегрирования

) Модель переводится в дискретный вид:

) Задается сигнал управления  и его производная

) Строится численная схема расчета с использованием схемы Эйлера.

Переходя к матричной форме:

Скрипт решения приведен ниже:

all;= 0; % Задание начальных условий= 1;

a21 = 2;= 3;= - 1;= - 2;= - 3;= - 4;= 0.02; % шаг

km = 100; % правая граница= zeros (2,km); % начальное обнуление Х(1,1) = 1; % начальные условия(2,1) = - 1;= [a11,a12; a21,a22];= [1,0; 0,1;];k=1: km %заполнение управляющих сигналов

kdt (k) = (k-1) *dt;(k) = k/200;(k) = 1;k=1: km-1(1) = b11*u1 (k) + b12*u2 (k);(2) = b21*u1 (k) + b22*u2 (k);(1,k+1) = dt*A (1,1) *x (1,k) + dt*A (1,2) *x (2,k) +x (1,k) + ub (1);(2,k+1) = dt*A (2,1) *x (1,k) + dt*A (2,2) *x (2,k) +x (2,k) + ub (2);

end= 1: 1: km;

%=========================Построение графиков('name','Управляющие воздействия u1 (k*dt),u2 (k*dt) '); %открытие окна для графика(kdt,u1 (k),kdt,u2 (k)); %рисуем график

set (gca,'YLim', [0 2]);('u1 (k*dt) ','u2 (k*dt) ');on; %сетка отображается('kdt','fontname','arial','fontsize',21,'color', [1 1 0]); % подпись Х('y','fontname','arial new','fontsize',23,'color', [0 1 0]); %подпись Y

figure ('name','Переходные процессы x1 (k*dt),x2 (k*dt) '); %открытие окна для графика(kdt,x (1,:),kdt,x (2,:)); %рисуем график

legend ('x1 (k*dt) ','x2 (k*dt) ');on; %сетка отображается('kdt','fontname','arial','fontsize',21,'color', [1 1 0]); % подпись Х('y','fontname','arial new','fontsize',23,'color', [0 1 0]); %подпись Y

Результаты работы приведены на Рис.3.1 и Рис.3.2 К отчету прилагается файл "Lab3_1. m" с реализацией этого задания.

Рис.3.1 - Управляющие воздействия

Рис.3.2 - Переходные процеессы

Сделаем идентификацию по схеме МНК.

Записываем исходную модель (в дискретной форме):

По условию измеряется , и переменная  не доступна к измерению. Она исключается из системы, для этого из первого уравнения выражается  и подставляется во второе:

И группируем переменные

) В качестве выходной переменной примем сигнал , переносим его в левую часть:

идентификация динамическая система сигнал

и окончательно (умножая правую часть на коэффициент ):

получим искомую линейно-регрессионную форму

,

где  - вектор параметров,

 - вектор регрессионных переменных.

) Формируется матрица  и вектор  для схемы МНК. Находится решение схемы МНК - оценку вектора параметров по уравнению:

Скрипт приведен ниже:

all;= 0; % Задание начальных условий= 1;

a21 = 2;= 3;= - 1;= - 2;= - 3;= - 4;= 0.02; % шаг

km = 100; % правая граница= zeros (2,km); % начальное обнуление Х(1,1) = 1; % начальные условия(2,1) = - 1;= [a11,a12; a21,a22];= [1,0; 0,1;];k=1: km %заполнение управляющих сигналов

kdt (k) = (k-1) *dt;(k) = k/200;(k) = 1;k=1: km-1(1) = b11*u1 (k) + b12*u2 (k);(2) = b21*u1 (k) + b22*u2 (k);(1,k+1) = dt*A (1,1) *x (1,k) + dt*A (1,2) *x (2,k) +x (1,k) + ub (1);(2,k+1) = dt*A (2,1) *x (1,k) + dt*A (2,2) *x (2,k) +x (2,k) + ub (2);

end= 1: 1: km;

%===========================================Метод МНК

%Формируем необходимые векторыk=1: km-1

dx (1,k) = (x (1,k+1) - x (1,k)) /dt;(1,km) = 0;(2,k) = (x (2,k+1) - x (2,k)) /dt;(2,km) = 0;(1,k) = (u1 (k+1) - u1 (k)) /dt;(1,km) = 0;(2,k) = (u2 (k+1) - u2 (k)) /dt;(2,km) = 0;= [dx (1,:) x (1,:) du (1,:) u1 du (2,:) u2];= vt';= vt*v;= v*x (1,:);

c1 = inv (F) *b;

% Матрица для проверки

c1 = [ (a11+a22) /a12 (a12*a21-a11*a22) /a12 b11/a12 (a22*b11-b21*a12) /a12 b12/a12 (a12*b22-b12*a22) /a12];

c1

Результаты работы приведены на Рис. 3.3

Рис.3.3 - Нахождение коэффициентов МНК

Сравнение

Сравнение параметров модели и результатов идентификации. Как видно из графика на Рис.3.4 параметры модели и результатов идентификации совпадают, что говорит о правильности решения.

Рис.3.4 - Сравнение параметров модели и результатов идентификации

Сравнение графиков модельного выходного сигнала и оценки выходного сигнала, восстановленной по схеме МНК

В качестве исходных данных для схемы МНК нам задан график (переходный процесс) выходного сигнала . Мы можем построить оценку этого сигнала , используя найденную оценку вектора параметров  по линейно-регрессионной формуле:

Сравнить модельный сигнал  и оценку сигнала  можно по формуле для суммы СКО для всех точек измерений:

При точных оценках вектора  СКО должно равняться нулю и графики  и  должны совпадать.

В построенной модели величина отклонения равна 0,02749, что подтверждает правильность решения.

Полный скрипт приведен ниже:

К отчету прилагается файл "Lab1_3. m" с реализацией этого задания

all;= 0; % Задание начальных условий

a12 = 1;= 2;= 3;= - 1;= - 2;= - 3;= - 4;= 0.02; % шаг

km = 100; % правая граница= zeros (2,km); % начальное обнуление Х(1,1) = 1; % начальные условия(2,1) = - 1;= [a11,a12; a21,a22];= [1,0; 0,1;];k=1: km %заполнение управляющих сигналов

kdt (k) = (k-1) *dt;(k) = k/200;(k) = 1;k=1: km-1(1) = b11*u1 (k) + b12*u2 (k);(2) = b21*u1 (k) + b22*u2 (k);(1,k+1) = dt*A (1,1) *x (1,k) + dt*A (1,2) *x (2,k) +x (1,k) + ub (1);(2,k+1) = dt*A (2,1) *x (1,k) + dt*A (2,2) *x (2,k) +x (2,k) + ub (2);

end= 1: 1: km;

%=========================Построение графиков('name','Управляющие воздействия u1 (k*dt),u2 (k*dt) '); %открытие окна для графика(kdt,u1 (k),kdt,u2 (k)); %рисуем график

set (gca,'YLim', [0 2]);('u1 (k*dt) ','u2 (k*dt) ');on; %сетка отображается('kdt','fontname','arial','fontsize',21,'color', [1 1 0]); % подпись Х('y','fontname','arial new','fontsize',23,'color', [0 1 0]); %подпись Y

figure ('name','Переходные процессы x1 (k*dt),x2 (k*dt) '); %открытие окна для графика(kdt,x (1,:),kdt,x (2,:)); %рисуем график

legend ('x1 (k*dt) ','x2 (k*dt) ');on; %сетка отображается('kdt','fontname','arial','fontsize',21,'color', [1 1 0]); % подпись Х('y','fontname','arial new','fontsize',23,'color', [0 1 0]); %подпись Y

%===========================================Метод МНК

%Формируем необходимые векторыk=1: km-1

dx (1,k) = (x (1,k+1) - x (1,k)) /dt;(1,km) = 0;(2,k) = (x (2,k+1) - x (2,k)) /dt;(2,km) = 0;(1,k) = (u1 (k+1) - u1 (k)) /dt;(1,km) = 0;(2,k) = (u2 (k+1) - u2 (k)) /dt;(2,km) = 0;= [dx (1,:) x (1,:) du (1,:) u1 du (2,:) u2];= vt';= vt*v;= v*x (1,:);

c1 = inv (F) *b;

% Сравнение параметров модели и результатов идентификации

c = [ (a11+a22) /a12 (a12*a21-a11*a22) /a12 b11/a12 (a22*b11-b21*a12) /a12 b12/a12 (a12*b22-b12*a22) /a12];

c= c;

% сравнение графиков модельного выходного сигнала и оценки

%выходного сигнала, восстановленной по схеме МНК

x11=c1*v;= 0;k=1: km= Otkl + (x (1,k) - x11 (k). ^2) / (x (1,k). ^2);

end