Построение математических моделей методом идентификации

Построение математических моделей методом идентификации

38

• Реферат
• Содержание
• Задание 1. Идентификация объектов методом наименьших квадратов
• Задание 2. Взаимосвязь различных форм моделей динамических систем
• Задание 3. Построение модели с распределенными параметрами
• Задание 4 численные процедуры оценивания параметров нелинейных регрессионных моделей
• Задание 5 разработка аналитических моделей объектов автоматизации. Линеаризация моделей
• Заключение

1.1 Постановка задачи

Математическая модель объекта имеет вид:

.(1.1)

Имеются статистические данные наблюдений:

Известна диагональная корреляционная матрица ошибок наблюдений

.

Требуется:

1. Оценить вектор коэффициентов модели В методом наименьших квадратов.

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

3. Оценить вектор коэффициентов модели методом максимального правдоподобия.

4. Построить доверительные интервалы для коэффициентов модели для доверительной вероятности 0,95.

5. Проверить значимость полученной статистической модели объекта.

6. Проанализировать полученные результаты.

Таблица 1.1

Вариант задания

1.2 Математическая постановка задачи

Пусть имеются транспонированные матрицы результатов наблюдений и след диагональной корреляционной матрицы ошибок D:

Поскольку результаты наблюдений суть случайные величины, получить «истинные» значения коэффициентов из модели

;(1.2)

нельзя. Вместо этого можно получить их оценки на основе исследований по таблице предварительно составленных наблюдений. Если речь идет о модели (1.2), то она принимает вид

,(1.3)

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

Из-за действия случайных возмущений предсказанное значение будет отличаться от результата измерения . Разности , , называют остатками.

Так как истинное значение вектора коэффициентов в и его оценка b различны,

, а .

Отсюда вектор остатков:

.

Оценки коэффициентов регрессии естественно искать так, чтобы обеспечить наименьшие возможные остатки, но остатки многочисленны, поэтому нужна некоторая суммарная характеристика, которая должна зависеть от различий между измеренными и предсказанными значениями выходных характеристик в каждом опыте. Такую функцию обычно называют функцией потерь, или функцией риска. Для различных целей и условий исследования она может иметь разный вид.

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

.(1.4)

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

Запишем сумму (1.4) в векторной форме. Пусть

обозначает N-мерный вектор столбец измеренных значений отклика, а

N-мерный вектор соответствующих им предсказанных значений, наконец,

вектор-столбец остатков. Как известно, скалярное произведение вектора на самого себя равно сумме квадратов его элементов, поэтому выражении (1.4) можно переписать в виде

.(1.5)

Метод, позволяющий оценивать регрессионные коэффициенты, выбирают так, чтобы минимизировать величину Q. Его называют обычно методом наименьших квадратов.

Пусть на основе данных таблицы исследования нужно найти такие оценки коэффициентов регрессии, которые минимизируют сумму Q, определенную в (1.4), минимум получим, приравняв производные по неизвестным оценкам к нулю. Но сначала подставим (1.2) в (1.4):

.(1.6)

После дифференцирования этого выражения по искомым оценкам и приравнивания нулю первых производных получаем систему уравнений:

(1.7)

.

Константы -2, входящие во все уравнения, не играют роли, ибо для равенства нулю произведения достаточно, чтобы равными нулю оказались соответствующие суммы. Поэтому полученная система сводится к виду:

(1.8)

Полученная система линейна относительно искомых оценок , а число уравнений в ней равно числу неизвестных коэффициентов k модели. Она называется системой нормальных уравнений.

Запись системы нормальных уравнений можно упростить, если положить

.(1.9)

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

(1.10)

В дальнейшем мы воспользуемся матричной записью. Обозначим - матрицу величин буквой G, - вектор оценок искомых коэффициентов буквой b, - вектор правой части системы буквой Z. Тогда

Важно заметить, что G - симметричная матрица, так как .

Вот матричная запись системы нормальных уравнений:

.(1.11)

Матрица G называется информационной матрицей. Ее можно представить через введенную в матрицу (1.2.2) матрицу регрессоров . Это утверждение проверяется непосредственно транспонированием F и перемножением . Точно так же проверяется, что . Значит, систему нормальных уравнений можно переписать так:

.(1.12)

Если ранг F равен k, то ранг то же равенk, так как из теории матриц известно, что произведение матриц есть положительно определенная матрица.

При этих условиях можно получить матрицу, обратную к информационной. Обозначим ее . Поскольку , умножение выражения (1.12) слева на матрицу приводит к решению системы нормальных уравнений:

.(1.13)

Таким образом, решив матричное уравнение (1.13), мы получаем в векторе b значение коэффициентов регрессии.

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

Для взвешенного метода наименьших квадратов матричное уравнение имеет вид:

.(1.14)

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

,(1.15)

где V(b) - дисперсия коэффициентов, - коэффициент Стьюдента. Коэффициент Стьюдента определяется по таблицам. Для нашего случая он равен 2,5. Дисперсия коэффициентов определяется по следующей формуле:

,(1.16)

где S - дисперсия оценки. Определяется следующим образом:

,(1.17)

Таким образом получим значения коэффициентов регрессии в виде

.

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

Метод наименьших квадратов

Расчет доверительных интервалов

Метод наименьших квадратов рис 1.1

МНК график остатков рис 1.2

2.1 Постановка задачи

Объект описан дифференциальным уравнением:

Требуется:

1. Записать модель объекта в пространстве состояний.

2. Записать модель объекта в форме передаточной функции.

3. Получить частотные характеристики объекта.

2.2 Математическая постановка задачи

Рассмотрим систему автоматического управления (САУ), описываемую линейным (линеаризованным) дифференциальным уравнением вида:

(2.1)

где u(t) - входной процесс, y(t) - выходной процесс, ai, bj - постоянные коэффициенты, n, m (n>m)- постоянные числа. В операторной форме выражение (2.1) может быть записано 38

.

Здесь D - оператор дифференцирования 38

. Отсюда преобразование “вход-выход” системы:

где W(D) называется операторной передаточной функции.

Один из способов моделирования систем заключается в представлении преобразования “вход-выход” в виде комплексной передаточной функции:

которая получается путем применения преобразования Лапласа к (2.2) ри начальных нулевых условиях. Здесь s-комплексная переменная. Связь между операторной (2.2) и комплексной (2.3) передаточными функциями можно записать в виде:

Комплексные числа, являющиеся корнями многочленаВ(s), называются нулями передаточной функции, а корни многочлена A(s) - полюсами.

Явный вид связи входа и выхода определяется выражением:

где w(t) - оригинал (т.е. полученный с помощью обратного преобразования Лапласа) комплексной передаточной функции W(s).

Динамические свойства систем характеризуют реакции на входные воздействия специального вида. В частности анализ выхода системы на единичный скачок и -функцию (дельта-функцию).

Пусть u(t) = 1(t), то есть на вход системы подается функция Хевисайда (единичный скачок), определяемая:

График функции Хевисайда приведен на рис. 2.1а:

а)

б)

Рис.2.1. Функции Хевисайда (а) и Дирака (б)

Реакция САУ на единичный скачек называется переходной функцией системы и обозначается h(t).

Если u(t) = (t), то есть на вход системы поступает функция Дирака (-функция, импульсная функция, рис. 2.1б) определяемая:

то реакция САУ называется импульсной переходной функцией системы и обозначается w(t). Таким образом оригинал комплексной передаточной функции можно измерить как реакцию систему на импульс.

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

Благодаря широкому применению при исследовании устойчивости динамических систем и проектировании регуляторов получили распространение частотные характеристики.

Пусть на вход системы с передаточной функцией W(s) подается гармонический сигнал u(t) = aucos(wt), t>0. В этих условиях справедлива следующая теорема:

Если звено является устойчивым, то установившаяся реакция y(t) на гармоническое воздействие является функцией той же частоты с амплитудой ay = au |W(iw)| и относительным сдвигом по фазе y = argW(iw).

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

y(t) = au |W(iw)| cos(w t + argW(iw)),(2.9)

где i - комплексная единица, - частотная характеристика.

При фиксированном значении w частотная характеристика является комплексным числом, и, следовательно, может быть представлена в виде:

где 38

- амплитудно-частотная характеристика (АЧХ); - фазово-частотная характеристика (ФЧХ); - вещественная частотная характеристика (ВЧХ); - мнимая частотная характеристика (МЧХ).

Геометрическое место точек W(iw) на комплексной плоскости при изменении w от w0 до от w1 (обычно w0 = 0, w1 = 38

), называется амплитудно-фазовой характеристикой (АФХ) или частотным годографом Найквиста.

Имеет широкое практическое значение диаграмма Боде (логарифмическая амплитудная характеристика, ЛАХ), которая определяется как L = 20 lg A(w), измеряется в децибелах и строится как функция от lg w .

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

Выполняем анализ с помощью пакета прикладных программ MatLab

Передаточная функция имеет вид:

Полюса передаточной функции:

-13.8796

-0.5602 + 1.4614i

-0.5602 - 1.4614i

Нули передаточной функции:

-35.7482

-0.2518

При помощи команды step(w) строим переходную функцию h(t) рис2.2

Рис.2.2. Переходная функция h(t)

При помощи команды impulse(w) строитьсяимпульсная переходная функция w(t) рис2.3:

Рис. 2.3. Импульсная переходная функцияw(t)

При помощи команды bode(w) строим Диаграмму Боде рис 2.4:

Рис.2.4. Диаграмма Боде

АФХ системы он жегодограф НайквистаW(i), строится при помощи команды nyquist(w):

Рис. 2.5. Годограф Найквиста

Текст расчетов в MatLab

>> w= TF([1,36,9],[1,15,18,34])

Transfer function:

s^2 + 36 s + 9

------------------------

s^3 + 15 s^2 + 18 s + 34

>>pole(w)

ans =

-13.8796

-0.5602 + 1.4614i

-0.5602 - 1.4614i

>>zero(w)

ans =

-35.7482

-0.2518

>>step(w)

>>impulse(w)

>>bode(w)

>>nyquist(w)

>>ltiview(w)

Рисунок 2.6 - Блок-схема программы исследования характеристик динамической системы

Плоская стенка первоначально прогрета равномерно до температуры 600С. В дальнейшем на внутренней поверхности стенки (х = 0) обеспечивается условие теплоизоляции (плотность теплового потока равна нулю), а с наружной поверхности (х = L) идет теплообмен с внешней средой, имеющей постоянную температуру Тср = -40 0С.. Изменение температуры в стенке осуществляется в результате процесса теплопроводности. Требуется получить зависимость от времени температуры на внутренней поверхности стенки. Толщина стенки L=20смб коэффициент теплоотдачи б=100Вт/м2К (неявная схема).

3.2 Математическая постановка задачи

Для неявной разностной схемы апроксимация уравнения теплопроводности будет иметь следующий вид:

(3.1)

Теперь удобно ввести и в левой и правой части сгруппировать все члены с индексрмj+1:

(3.2)

Известная функция u(i,j) связана функциональной зависимостью с тремя неизвестными функциями на последующем слое. Для заданного временного слоя такие уравнения надо записать для всех внутренних точек i=2….n-1 и дополнить их граничными условиями i=1,n. Тогда получим систему с трех диагональной матриц неизвестных членов, которая решается с помощью метода прогонки.

Так как в задаче заданна температура окружающей среды, то имеем граничные условия 3-го рода. Запишем уравнения для граничных и внутренних точек одного временного слоя.

Для i=1, имеем:

(3.3)

(3.4)

Обозначим , , получим уравнение для i=1:

(3.5)

Для i=2…n-1, имеем:

(3.6)

(3.7)

(3.8)

Для i=n, имеем:

(3.9)

(3.10)

Обозначим , , получим уравнение для n=1:

(3.11)

Рассмотрим метод прогонки для решения системы, сотоящей из уравнений (3.5), (3.8), (3.11). Из уравнения (3.5) выразим U1,j+1:

(3.12)

, (3.13)

Запишем уравнение для i=2 (3.8) и подставим в него выражение (3.12) и выразим U2,j+1:

(3.14)

(3.15)

(3.16)

Продолжим этот процесс до i=n-1, получим:

(3.17)

Коэффициенты fi и giизвестны из граничных условий на первой границе f1 =U1,j и gi=0, их называют прогоночными коэффициентами, и мы можем их найти по возрастающей рекурсии вплоть до i=n-1. Можем записать:

(3.18)

Подставим выражение (3.18) в уравнение для второй границы (3.11) и выразим Un,j+1:

(3.19)

(3.20)

Теперь Un,j+1нам известна. Используя полученное выражение, определим Un-1,j+1из выражения (3.18) для всех i, включая i=1. Таким образомUш,j+1 мы определяем по обратной рекурсии от i+1 к i, в то время как fi и giопределяли по прямой рекурсии отi-1 к i. Такой метод расчета дает наименьшую погрешность.

3.3 Описание входных и выходных данных

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

3.4 Текст программы в SciLab

//Заданные значения

t0=873;

L=0.2;

alfa=100;

lamda=58.7;

c=500;

p=7900;

tos=233;

//Количество шагов и шаг по длинe стержня

n=20;

dx=(L-0)/n;

//Шаг по времени, время,кол-во шагов и сетка по времени

dt=100;

time=70000;

m=time/dt;

T=0:dt:time;

//Начальное распределение температур по стержню

for i=1:n+1

U(i,1)=t0;

end

//Температура по длине стержня во времени по неявной схеме

l=dt/dx^2;

k=lamda/c/p;

b=alfa/c/p;

c=dt/dx;

g(1)=0;

for j=2:m+1

f(1)=U(1,j-1);

for i=2:n+1

f(i)=(k*l*f(i-1)+U(i,j-1))/(1+2*k*l-k*l*g(i-1));

g(i)=k*l/(1+2*k*l-k*l*g(i-1));

end

U(n+1,j)=(2*b*c*tos+2*k*l*f(n)+U(n+1,j-1))/(1+2*k*l+2*b*c-2*k*l*g(n));

for i=n:-1:2

U(i,j)=g(i)*U(i+1,j)+f(i);

end

U(1,j)=(U(1,j-1)+2*k*l*U(2,j))/(1+2*k*l);

end

//Температура изолированного конца стержня

scf(1);

plot(T,U(n+1,:),'k-');

xgrid(1);

3.5 Результаты работы программы

3.5 Блок-схема программы

Рисунок 3 - Блок-схема программы для расчета температурного режима плоской стенки

4.1 Постановка задачи

Пусть некоторый процесс описывается внутренне нелинейной моделью, т.е. такой, которая не может быть представлена в линейном виде относительно параметров. Будем искать уравнение модели в виде:

(4.1)

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

4.2 Математическая постановка задачи

В методе наискорейшего спуска направление движения при поиске минимума F(x) задается вектором антиградиента ? gradF(p1,p2) функции F(p1,p2) в рассматриваемой точке, т.е:

( (4.2)

Функция F(p1,p2) расчитывается как квадрат ошибки расчетных и фактических значений отклика:

(4.3)

Составляющие вектора градиента в точке xk определяются значениями частных производных первого порядка функции F(x) по соответствующим переменным, вычисленным в точке (p1k , p2k):

(4.4)

Направление вектора градиента совпадает с направлением наискорейшего возрастания функции F(p1,p2). Вектор ? gradF(p1,p2) называется антиградиентом, его направление совпадает с направлением наискорейшего убывания функции. Предполагается, что компоненты градиента могут быть записаны в аналитическом виде или с достаточно высокой точностью вычислены при помощи численных методов.

Выбор величины шага л осуществляется путем решения задачи минимизации F(p1,p2) в направленииgradF(p1,p2) с помощью одного из методов одномерного поиска.

Преимущество метода в том, что при переходе от шага к шагу обеспечивается выполнение неравенстваF(p1 k+1,p2k+1)<= F(p1 k,p2k), т.е. значение функции цели улучшается.

Скорость сходимости метода при решении некоторых задач является недопустимо низкой.

Это связано с тем, что изменения переменных непосредственно зависят от величины градиента, которая стремится к нулю в окрестности точки минимума. Поэтому метод наискорейшего спуска эффективен при поиске на значительных расстояниях от точки минимума x* и плохо работает в окрестности этой точки.

4.3 Описание входных и выходных переменных

В качестве входных переменных имеем массив независимых значений Х и фактические значения отклика YF, точность вичислений eps.

Выходные величины - конечне значения p1 и p2, при которых достигается минимум функции F(p1,p2), количество итераций, что характеризует скорость сходимости метода и значение функции F(p1,p2), а также график сравнения фактического и расчетного значения отклика и таблицу с фактическим и расчетным значениями отклика.

4.4 Текст программы в SciLab

//Входные данные

X1=[1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6];

YF=[20 30 70 85 100 90 100 108 120 110 124]; //Y фактическое

//Расчетная функции в виде у=р1-р2expХ1)

function YR=raschet(P)

YR=P(1)-P(2)*exp(-X1);

endfunction

//Минимизируемая функция как результат

//ошибки расчетных и фактических значений

function F=oshibka(P)

F=sum((YF-raschet(P))^2);

endfunction

eps=0.00001; //точность вычислений

P=[10 0.9]; //начальное приближение

lamda=1; //параметр,характеризующий длину шага

F0=oshibka(P*2);

F=oshibka(P);

i=0;

while abs(F0-F)>eps

i=i+1;

P0=P;

GR=numdiff(oshibka,P0);

P=P0-lamda*GR;

F0=oshibka(P0);

F=oshibka(P);

if F0<F

lamda=lamda/2;

end

end

Значения р1 и р2, при которых ошибка расчетов минимальна'

p1=P(1)

p2=P(2)

Скорость схождения(количество циклов)'

i

Y=raschet(P);

'Таблица результатов расчета'

[X1;YF;Y;(YF-Y)^2]'

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

sigma=F

scf(1);

plot(X1,YF,'k.',X1,Y,'k-');

xgrid(1);

4.5 Результаты работы программы

Значения р1 и р2, при которых ошибка расчетов минимальна

p1 = 110.90905

p2 = 282.44884

Скорость схождения(количество циклов)

i = 502.

Таблица результатов расчета

XYFYR (YF-YR)^2

1. 20. 7.0019302 168.94982

1.5 30. 47.886197 319.91604

2. 70. 72.683758 7.2025575

2.5 85. 87.724239 7.4214795

3. 100. 96.846752 9.9429717

3.5 90. 102.37984 153.26034

4. 100. 105.73582 32.899642

4.5 108. 107.77133 0.0522905

5. 120. 109.00593 120.86965

5.5 110. 109.75475 0.0601485

6. 124. 110.20893 190.19358

Среднее квадратическое отклонениеsigma = 1010.7685

4.6 Блок-схема программы

Разработайте математическую модель заданной динамической системы. Построить аналитическую модель емкости с одним притоком и двумя стоками (например двухручьевая МНЛЗ) и провести ее линеаризацию. Математическую модель заданного объекта представить в виде программы на ЭВМ, рассчитать кривые исходной и линеаризованной модели отклика модели на ступенчатое и импульсное воздействие. Предложить методику идентификации параметров модели. Какие эксперименты будет необходимо провести на реальном объекте? Оцените необходимый объем и форму представления результатов.

5.2 Математическая постановка задачи

В данном объекте управления выходной величиной является уровень металла в промежуточном ковше h, а входной величиной - разность между притоком и стоками металла ?G=Gпр-Gст1-Gст2, причем возмущения могут возникать за счет изменения Gпр,Gст1 и Gст2.

Изменение уровня характеризуется следующим дифференциальным уравнением:

, (5.1)

где с - плотность металла, F - площадьзеркала металла в промковше.

Известно, что сток жидкости через отверстие пропорционален корню квадратному из высоты этой жидкости на отверстием:

(5.2)

где бпр,бст1,бст2 - коэффициенты расхода на притоке и стоках; fпр,fст1,fст2 - проходные сечения в стопорных парах или в шиберных затворах сталеразливочного и промежуточного ковшей; H(t) - текущее значение уровня металла в сталеразливочном ковше.

С использование зависимостей (5.2) уравнение (5.1) приобретает вид нелинейного дифференциального уравнения первого порядка:

(5.3)

Уравнение можно несколько упростить, подвергнув линеаризации Gпр,Gст1 и Gст2, использовав разложение в ряд Тейлора Gпр в окрестностях точки Н0 и Gст1 и Gст2 в окрестностях точки h0 и отбросив члены ряда, содержащие величины второго и более порядков малости.

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

(5.4)

где H и h - уровни металла в сталеразливочном и промежуточном ковшах;

Т - постоянная времени объекта;

k1,k21,k22,k3 - коэффициенты передачи объекта по различным каналам возмущений.

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

(5.5)

Таким образом правая часть уравнения (5.4) учитывает четыре возможных возмущения: по каналам регулирующих органов на притоке и стоках и каналу изменения уровня металла в сталеразливочном ковше.

Решая уравнения (5.3) и (5.4) с помощью конечных разностей, получим следующие выражения для нелинеаризованной и линеаризованной моделей:

(5.6)

(5.7)

5.3 Результаты расчетов модели промежуточного ковша двухручьевой МНЛЗ

а) нелинеаризованная модель (формула (5.6)):

- установившийся режим

- единичное воздействие (функция Хевисайда)

- импульсное воздействие (функция Дирака)

б) линеаризованная модель (формула (5.7)):

- единичное воздействие (функция Хевисайда)

- импульсное воздействие (функция Дирака)

Графики построены для возмущения по каналу уровня металла в сталеразливочном ковше.