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

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

Министерство образования и науки, молодежи и спорта Украины

Государственное высшее учебное заведение

Приазовский государственный технический университет

Факультет информационных технологий

Кафедра автоматизации и компьютерных технологий

ПОЯСНИТЕЛЬНАЯ ЗАПИСКА

по курсовой работе по дисциплине  “Идентификация и моделирование  объектов автоматизации”

Выполнил:  студентка группы МА-10  Адамова О.С.

Мариуполь 2013 г.

Реферат

Пояснительная записка содержит 42 страницы, 12 рисунков, 3 таблицы.

Целью работы является получение численного решения для конкретной задачи. Это решение должно быть получено в результате выполнения соответствующей программы на ЭВМ, написанной на языке высокого уровня, составленной самим обучающимся.

При выполнении работы необходимо предварительно ознакомиться с соответствующим методом и его алгоритмом. Выполнение программы на ЭВМ и анализ получения результатов составляют заключительную стадию работы. Реализация программы должна быть показана на контрольном примере с выдачей выходных документов на принтер.

Введение

Цель работы − получение практических навыков в построении математических моделей технических объектов, написании программ для решении задач моделирования с использованием языка программирования С/С++ и математических пакетов MathCad или MatLab, изучение теоретических основ и особенностей выполнения параметрической идентификации различных моделей, реализации алгоритмов линейного и нелинейного регрессионного анализа, планирования эксперимента.

Задачи курсовой работы включают:

− получение студентами навыков самостоятельной работы;

− освоение технологии разработки и отладки программ, реализующих модели технических объектов;

− более качественное изучение нормативных материалов - государственных стандартов и технических условий;

− более полное изучение базовых средств языков программирования и получение навыков постановки и решения различных задач с помощью ПЭВМ;

− изучение и использование сред численного моделирования и статистического анализа (MatLab, StatGraph и т.п.).

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

В ходе выполнения курсовой работы необходимо разработать программы на Scilab, Matlab, C++, которая позволит:

− оценить построенную математическую модель;

− найти выходные параметры, описывающие математическую модель, и построить зависимости между входными и выходными характеристиками объекта;

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

Входные и выходные данные

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

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

1. Идентификация объектов методом наименьших квадратов

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

Матрица X

Матрица Y

Для линейных уравнений вида:

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

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

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

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

. Найти скорректированный коэффициент множественной детерминации. Сравнить его с нескорректированным (общим) коэффициентом детерминации.

4. С помощью F-критерия Фишера оценить статистическую надежность уравнения регрессии и коэффициента детерминации

. С помощью t-критерия Стьюдента оценить статистическую значимость коэффициентов чистой регрессии.

. Составить уравнение линейной парной регрессии, оставив лишь один значащий фактор.

Для наших данных система нормальных уравнений имеет вид:

Уравнение множественной регрессии:

Оценка значимости уравнения регрессии в целом производится не основе  F-критерия Фишера, которому предшествует дисперсионный анализ. Согласно основной идее дисперсионного анализа, общая сумма квадратов отклонений переменной  от среднего значения  раскладывается на две части - объясненную и необъясненную:

 - общая сумма квадратов отклонений;

 - сумма квадратов отклонений, объясненная регрессией (факторная сумма квадратов отклонений);

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

Сопоставляя факторную и остаточную дисперсии в расчете на одну степень свободы, получим величину F-критерия Фишера:

Фактическое значение F-критерия Фишера сравнивается с табличным значением  при уровне значимости  и степенях свободы   и  При этом, если фактическое значение F-критерия Фишера больше табличного, то уравнение признается статистически значимым:

Качество модели, исходя из относительных отклонений по каждому наблюдению, признается хорошим, т.к. средняя ошибка аппроксимации не превышает 10 %.

Средние коэффициенты эластичности для линейной регрессии:

Средние коэффициенты эластичности показывают, на сколько процентов в среднем изменится результат при изменении соответствующего фактора на 1 %. Таким образом, подтверждается большее влияние на результат  фактора  чем факторов

Показателем интенсивности связи служит значение коэффициента корреляции. Считается, если он равен 1, то взаимозависимость признаков является строгой (полной); если его значение находится в интервале от 1 до 0,8, то это свидетельствует о сильной их взаимозависимости; если в интервале от 0,7 до 0,3 - об умеренной (не ярко выраженной) взаимозависимости, а если же оно лежит в интервале от 0,2 до 0,0, то мы имеем дело со слабой или нулевой взаимозависимостью.

Коэффициенты парной корреляции:

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

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

− определитель матрицы парных коэффициентов корреляции

− определитель матрицы межфакторной корреляции

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

Эта доля составляет 99,9 % и указывает на высокую степень обусловленности вариации результата вариацией факторов (тесную связь факторов с результатом).

Скорректированный коэффициент множественной детерминации определяет тесноту связи с учетом степеней свободы общей и остаточной дисперсий:

Оба коэффициента указывают на высокую детерминированность результата  в модели факторами

Оценка статистической значимости параметров регрессии проводится по  t-критерию Стьюдента:

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

2. Исследование разомкнутой линейной системы

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

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

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

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

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

. Определить нули и полюса передаточной функции.

. Определить коэффициент усиления системы в установившемся режиме и полосу пропускания системы.

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

. Построить процесс на выходе системы при произвольном входном сигнале.

. Использовать модуль LTI-Viewer для построения различных характеристик.

Для описания линейных систем могут применяться несколько способов:

− дифференциальные уравнения;

− модели в пространстве состояний;

− передаточные функции;

− модели вида «нули-полюса».

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

Модель объекта в форме передаточной функции:

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

all;;

% Ввод передаточной функции %

num=[0 12 6 5]=[1 3 2 -4]=tf(num,den)

% Построение модели объекта в пространстве состояний %_ss=ss(w)

% Нахождение нулей и полюсов передаточной функции %

z=zero(w)=pole(w)

% Нахождение коэффициента усиления системы %

% в установившемся режиме %=dcgain(w)

% Определение полосы пропускания системы %=bandwidth(w)

% Построение модели системы в форме "нули-полюса" %_zpk=zpk(w)

% Расположение нулей и полюсов системы на графике %

pzmap(w);grid;-dmeta;

% Построение переходной функции %(w);grid;-dmeta;

% Построение импульсной переходной функции %

impulse(w);grid;-dmeta;

% Создание массива частот для построения %

% амплитудно-частотной характеристики %=logspace(-4,4,500);

r=freqresp(w,freq);=r(:);(freq,abs(r));grid;

print -dmeta;

% Создание массива частот для построения %

% фазо-частотной характеристики %=logspace(-4,4,500);

r=freqresp(w,freq);=r(:);=angle(r)*180/pi;(freq,phi);grid;-dmeta;

% Диаграмма Боде %(w);grid;

print -dmeta;

% Частотный годограф Найквиста %

nyquist(w);grid;-dmeta;

% Сигнал, имитирующий прямоугольные импульсы %

% единичной амплитуды %

% (период - 4 секунды, количество - 5 импульсов) %

[u,t]=gensig('square',4);

lsim(w,u,t);grid;

print -dmeta;

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

= 0 12 6 5= 1 3 2 -4= 12 s^2 + 6 s + 5

--------------------^3 + 3 s^2 + 2 s - 4time transfer function._ss = = x2 x3-3 -1 22 0 00 1 0= 400= x2 x33 0.75 0.625= 0time state-space model.=

.2500 + 0.5951i

.2500 - 0.5951i=

.8982 + 1.1917i

.8982 - 1.1917i

.7963 + 0.0000i=

.2500=

.3819_zpk =

(s^2 + 0.5s + 0.4167)

(s-0.7963) (s^2 + 3.796s + 5.023)time zero/pole/gain model.

Рис. 2.1. Расположение нулей и полюсов системы на графике

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

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

Рис. 2.4. Амплитудно-частотная характеристика

Рис. 2.5. Фазо-частотная характеристика

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

Рис. 2.7. Частотный годограф Найквиста

Рис. 2.8. Сигнал, имитирующий прямоугольные импульсы единичной амплитуды (период − 4 секунды,  количество − 5 импульсов)

3. Построение модели с распределенными параметрами

.1 Исходные данные

Рассмотрим стержень из теплопроводящего материала с коэффициентом теплопроводности k. Предположим, что температура на концах стержня задана, а боковая поверхность стержня теплоизолирована. Пусть ось x направлена вдоль оси стержня, а его концы расположены в точках x=0 и x=L.

.2 Общие теоретические сведения

Задача сводится к определению зависимости от времени t температуры u в точках стержня, то есть функции двух переменных u(x,t). Функция u(x,t) должна удовлетворять уравнению теплопроводности

 (0<x<L),(3.1)

начальному условию

u(x,0)=f(x), (0<x<L),(3.2)

и условиям на концах стержня

u(0,t)=j₁(t), u(L,t)=j₂(t), (tV0).(3.3)

Значения u(0,0) и u(L,0), полученные из (2) и (3), должны совпадать. Это будет если jj₁(0)=f(0), jj₂(0)=f(L).

Следует отметить, что путем замены переменных t^ў=a²t уравнение (1) можно преобразовать к виду:

(3.4)

Это означает, что решение задачи (3.1÷3.3) путем замены переменных сводится к решению задачи (3.4). Далее будем полагать а=1.

Построим на плоскости (x,t) сетку с шагом h по переменной x  (x_i=(i-1)h; i=1, …, n+1; h=L/n) и с шагом tt по переменной t (t_j=(j-1)tt). Обозначим u_ij=u(x_i,t_j). Производные в уравнении (1) аппроксимируем следующим образом:

(3.5)

(3.6)

Подставляя (3.5) и (3.6) в (3.1) при a=1, получим разностное уравнение:

(3.7)

В соответствии с (3.2) и (3.3) значения

u_i0=f(x_i), u_0j=j₁(t_j), u_nj=j₂(t_j)(3.8)

являются известными. Тогда, подставляя в (3.7) j=0, получим систему n-1 линейных уравнений, решив которую можно определить u_i1 (i=1, …, n-1). При этом, поскольку u₀₁=jj₁(t₁), …, u_n1=jj₂(t₁), известными оказываются все значения временного слоя j=1, (t=t₁). Затем, подставляя в (3.7) j=2, решаем систему уравнений относительно u_i2 и т.д. для всех j=2, …, m.

Из (3.7) следует, что в каждое i-тое уравнение (i=1, …, n-1) с ненулевыми коэффициентами входят только три неизвестных (u_i-1,j, u_ij, u_i+1,j). Величина u_i,j-1 к этому моменту является известной и потому отнесена в правую часть уравнения. Таким образом, матрица системы уравнений является трехдиагональной и эту систему можно решить методом прогонки. Для этого представим ее в стандартном виде:

(3.9)

Для данной задачи x_i=u_ij, a_i=l, g_i=l, b_i=1-2l, b₀=1, g₀=0, j₀=u_0j=j₁(t_j), j_n=u_nj=j₂(t_j), j_i=-u_i,j-1 (i=1, …, n-1).

Пусть на j-том шаге заданными являются параметры  u_i,j-1 (i=1, …, n-1), u_0j, u_nj, ll. Все неизвестные значения u_ij можно разместить в массиве x_i (x_i=u_ij; i=0, …, n). Ищем связь x_i-1 с x_i в виде рекуррентного соотношения:

x_i-1=c_i-1x_i+n_i-1, i=1, …, n.(3.10)

Подставляя (3.10) в (3.7), получаем:

lc_i-1x_i-(1+2l)x_i+lx_i+1 = -u_i,j-1-ln_i-1.(3.11)

Отсюда

(3.12)

Сравнивая (3.12) с (3.10), находим рекуррентные соотношения:

,

      (3.13)

cc₀=0, nn₀=u_0j.

Таким образом, алгоритм определения значений u_ij по известным u_i,j-1 состоит из двух этапов: прямого хода прогонки по формулам (3.13) при (i=1, …, n-1) и обратного хода прогонки по формуле (3.10) при (i=n, …, 2).

а)  б)

Рис. 3.1. Шаблоны неявной (а) и явной (б) разностных схем

Необходимо отметить, что разностное уравнение (7) связывает одно известное значение U_i,j-1 (из предыдущего j-1 временного слоя) и три неизвестных (U_i,j, U_i-1,j, U_i+1,j). Поэтому найти значения U_i,j (i=1, …, n-1) можно только все сразу путем решения системы уравнений. Такая схема связи переменных в разностном уравнении называется неявной. Шаблон неявной разностной схемы представлен на рис. 3.1 (а).

Наряду с неявной возможна организация явной разностной схемы. Для этого вместо выражения (5) для первой разностной производной по времени используют формулу:

(3.14)

Тогда разностное уравнение запишется в виде:

(3.15)

В этом случае связываются три неизвестные значения, относящиеся к предыдущему временному слою (здесь j-тому) и только одно неизвестное U_i,j+1. Шаблон явной разностной схемы представлен на рис. 3.1 (б).

При использовании этой схемы неизвестные параметры определяются путем последовательного применения формулы (2.14) при i=1, …, n-1. Поскольку при этом не надо решать системы уравнений, то процесс определения параметров одного временного слоя требует меньших затрат времени, чем в случае неявной схемы.

Однако, неявная схема устойчива (ошибка не возрастает от шага к шагу) при любых значениях λ=τ/h². Явная схема является устойчивой только при λ <1/2. В противном случае развивается экспоненциальный рост погрешности так, что обычно происходит аварийная остановка ЭВМ по переполнению порядка. Поэтому при использовании явной схемы вычисления приходится вести с очень малым шагом по времени.

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

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

j₁(t)Tf(0), j₂(t)Tf(L).

.3 Решение задачи для исходных данных

Решить смешанную задачу для уравнения теплопроводности с начальными u(x,0)=f(x) и граничными условиями u(0,t)=f(0), u(1,t)=b, L=1.

Табл. 3.1. Вариант задания

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

//---------------------------------------------------------------------------

#include <vcl.h>

#pragma hdrstop

#include "Unit1.h"

//---------------------------------------------------------------------------

#pragma package(smart_init)

#pragma resource "*.dfm"*Form1;

//---------------------------------------------------------------------------

__fastcall TForm1::TForm1(TComponent* Owner)

: TForm(Owner)

{

}

//---------------------------------------------------------------------------

#define a 14

#define b 16

#define c 0.25

#define d 0.5

#define N 20

#define L 1

#define tau 60

#define ro 7800

#define lam 46

#define c_f 460

#define t_max 3000T[N+1],alpha[N+1],beta[N+1],time1=0,h;A[N+1],C[N+1],B[N+1],F[N+1];h_iav,tau_iav,a_iav,T_iav[N+1],TT_iav[N+1],time_iav=0;funkf(float x)

{m;(x>=0)

{(x<c)

{=a-(x/c)*a;

}(x<d)

{=0;

}

else

if(x<=1)

m=((1-x)*b)/(1-d);

}(m);

}start_heat_field(float h)

{(int i=1;i<=N;i++)[i]=funkf(h*(i-1));

}neiavnii_metod_start(){=(float)L/(float)(N);_heat_field(h);(int i=1;i<=N;i++)->Series1->AddXY(h*i,T[i]);

}neiavnii_metod_calc(){=time1+tau;[1]=0;[1]=funkf(0);(int i=2;i<=N-1;i++){[i]=C[i]=lam/(h*h);[i]=(2*lam)/(h*h)+(ro*c_f)/tau ;[i]=-((ro*c_f)/tau)*T[i];[i]=A[i]/(B[i]-C[i]*alpha[i-1]) ;[i]=(C[i]*beta[i-1]-F[i])/(B[i]-C[i]*alpha[i-1]);

}[N]=funkf(L-0.0000001);(int i=N-1;i>=2;i--){[i]=(A[i]/(B[i]-C[i]*alpha[i]))*T[i+1]+(C[i]*beta[i-1]-F[i])/(B[i]-C[i]*alpha[i-1]);->Series1->Clear();[0]=funkf(0.0000001);(int i=1;i<=N;i++)->Series1->AddXY(h*i,T[i]);

}

}iavnii_metod_start(){_iav=lam/(ro*c);_iav=(h*h)/(4*a_iav);(int i=1;i<=N;i++)_iav[i]=funkf(h*(i-1));(int i=1;i<=N;i++)->Series2->AddXY(h*i,T_iav[i]);_iav[1]=funkf(0.0000001);_iav[N]=funkf(L-0.0000001);

}iavnii_metod_calc(){_iav=time_iav+tau_iav;(int i=1;i<=N+1;i++)_iav[i]=T_iav[i];(int i=2;i<=N-1;i++)_iav[i]=TT_iav[i]+((lam*tau_iav)/(ro*c))*((TT_iav[i+1]-2*TT_iav[i]+TT_iav[i-1])/(h*h));->Series2->Clear();(int i=1;i<=N;i++)->Series2->AddXY(h*i,T_iav[i]);

}__fastcall TForm1::FormCreate(TObject *Sender)

{_metod_start();_metod_start();

}

//---------------------------------------------------------------------------__fastcall TForm1::Button1Click(TObject *Sender)

{=0;_metod_start();_metod_start();(time1<t_max){_metod_calc();_metod_calc();

}->Caption="Time= "+FloatToStr(time1);

}

//---------------------------------------------------------------------------__fastcall TForm1::Button2Click(TObject *Sender)

{=0;_metod_start();_metod_start();(time1<TrackBar1->Position*tau){_metod_calc();_metod_calc();

}->Caption="Time= "+FloatToStr(time1);

}

//---------------------------------------------------------------------------

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

Вывод:

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

корреляция модель итерационный функция

4. Численные процедуры оценивания параметров нелинейных регрессионных моделей

Метод Ньютона, алгоритм Ньютона (также известный как метод касательных) − это итерационный численный метод нахождения корня (нуля) заданной функции. Поиск решения осуществляется путем построения последовательных приближений и основан на принципах простой итерации. Метод обладает квадратичной сходимостью. Улучшением метода является метод хорд и касательных. Также метод Ньютона может быть использован для решения задач оптимизации, в которых требуется определить нуль первой производной либо градиента в случае многомерного пространства.

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

Пусть  − определена на отрезке  и дифференцируемая <#"723183.files/image093.gif"> 

где  − угол наклона касательной в точке .

Следовательно, искомое выражение для  имеет вид:

Итерационный процесс начинается с некоего начального приближения  (чем ближе к корню, тем лучше, но если предположения о его нахождении отсутствуют, методом проб и ошибок можно сузить область возможных значений, применив теорему о промежуточных значениях <#"723183.files/image098.gif">.

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

Недостатки:

. Если начальное приближение недостаточно близко к решению, то метод может не сойтись.

. Если производная <#"723183.files/image102.gif">

Вывод:

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

Заключение

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

Перечень ссылок

1. Автоматизированные системы управления технологическими процессами. Идентификация и оптимальное управление / Под ред. Салнги В.И. − Харьков: «Вища школа», 1976 г. − 180 с.

. Брикман М.С., Кристинков Д.С. Аналитическая идентификация управляемых систем. − Рига: «Зинатне», 1974 г. − 206 с.

. Гельфандбейн Я.А. Методы кибернетической диагностики динамических систем. Идентификация функционирующих систем математическими моделями. − Рига: «Зинатне», 1967 г. − 542 с.

. Дрейпер Н., Смит Г. Прикладной регрессионный анализ. − М.: «Статистика», 1973 г. − 391 с.

. Демиденко Е.З. Линейная и нелинейная регрессии. - М.: «Финансы и статистика», 1981 г.