Идентификация и моделирование объектов автоматизации

print -dmeta;

 

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

 

num =

 

     0    12     6     5

 

den =

 

     1     3     2    -4

 

w =

 

    12 s^2 + 6 s + 5

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

  s^3 + 3 s^2 + 2 s - 4

 

Continuous-time transfer function.

 

w_ss =

 

  a =

       x1  x2  x3

   x1  -3  -1   2

   x2   2   0   0

   x3   0   1   0

 

  b =

       u1

   x1   4

   x2   0

   x3   0

 

  c =

          x1     x2     x3

   y1      3   0.75  0.625

 

  d =

       u1

   y1   0

 

Continuous-time state-space model.

 

z =

 

  -0.2500 + 0.5951i

  -0.2500 - 0.5951i

 

p =

 

  -1.8982 + 1.1917i

  -1.8982 - 1.1917i

   0.7963 + 0.0000i

 

k =

 

   -1.2500

 

b =

 

    0.3819

 

w_zpk =

 

      12 (s^2 + 0.5s + 0.4167)

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

  (s-0.7963) (s^2 + 3.796s + 5.023)

 

Continuous-time zero/pole/gain model.

 

 

 

 

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

 

 

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

 

 

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

 

 

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

 

 

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

 

 

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

 

 

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

 

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

 

 

 

 

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

 

3.1 Исходные данные

 

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

 

3.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)=jj1(t),   u(L,t)=jj2(t),   (tVł0).

(3.3)

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

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

(3.4)

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

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

	(3.5)
	(3.6)

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

(3.7)

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

ui0=f(xi),   u0j=jj1(tj),  unj=jj2(tj)

(3.8)

являются известными. Тогда, подставляя в (3.7)  j=0,   получим систему   n-1  линейных уравнений, решив которую можно определить  ui1   (i=1, …, n-1).   При этом, поскольку  u01=jj1(t1), …, un1=jj2(t1),   известными оказываются все значения временного слоя  j=1,  (t=t1).   Затем, подставляя в (3.7)  j=2,   решаем систему уравнений относительно  ui2   и т.д. для всех  j=2, …, m.

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

    (3.9)

Для данной задачи   xi=uij, aai=ll,   ggi=ll, bbi=1-2ll,  bb0=1, gg0=0,   jj0=u0j=jj1(tj), jjn=unj=jj2(tj),   jji=-ui,j-1   (i=1, …, n-1).

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

xi-1=cci-1xi+nni-1, i=1, …, n.

(3.10)

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

llcci-1xi-(1+2ll)xi+llxi+1 = -ui,j-1-llnni-1.

(3.11)

Отсюда

(3.12)

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

,

(3.13)

cc0=0,  nn0=u0j.

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

 

а)

б)

 

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

 

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

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

(3.14)

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

(3.15)

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

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

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

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

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

jj1(t)Tşf(0),  jj2(t)Tşf(L).

 

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

 

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

 

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

 

№	а	b	c	d	f(x)
1	1,1	3	0,05

 

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

 

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

 

#include <vcl.h>

#pragma hdrstop

 

#include "Unit1.h"

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

#pragma package(smart_init)

#pragma resource "*.dfm"

TForm1 *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 3000

 

float T[N+1],alpha[N+1],beta[N+1],time1=0,h;

float A[N+1],C[N+1],B[N+1],F[N+1];

float h_iav,tau_iav,a_iav,T_iav[N+1],TT_iav[N+1],time_iav=0;

 

float funkf(float x)

{

float m;

if(x>=0)

{

if(x<c)

{

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

}

else

if(x<d)

{

m=0;

}

                        else

                                if(x<=1)

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

        }

return(m);

}

 

void start_heat_field(float h)

{

        for(int i=1;i<=N;i++)

                T[i]=funkf(h*(i-1));

}

 

void neiavnii_metod_start(){

h=(float)L/(float)(N);

start_heat_field(h);

for(int i=1;i<=N;i++)

        Form1->Series1->AddXY(h*i,T[i]);

 

}

 

void neiavnii_metod_calc(){

time1=time1+tau;

alpha[1]=0;

beta[1]=funkf(0);

 

 

for (int i=2;i<=N-1;i++){

        A[i]=C[i]=lam/(h*h);

        B[i]=(2*lam)/(h*h)+(ro*c_f)/tau ;

        F[i]=-((ro*c_f)/tau)*T[i];

        alpha[i]=A[i]/(B[i]-C[i]*alpha[i-1]) ;

        beta[i]=(C[i]*beta[i-1]-F[i])/(B[i]-C[i]*alpha[i-1]);

}

T[N]=funkf(L-0.0000001);

for (int i=N-1;i>=2;i--){

        T[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]);

Form1->Series1->Clear();

 

T[0]=funkf(0.0000001);

for(int i=1;i<=N;i++)

        Form1->Series1->AddXY(h*i,T[i]);

 

}

}

 

void iavnii_metod_start(){

a_iav=lam/(ro*c);

 

tau_iav=(h*h)/(4*a_iav);

 

for(int i=1;i<=N;i++)

        T_iav[i]=funkf(h*(i-1));

 

for(int i=1;i<=N;i++)

        Form1->Series2->AddXY(h*i,T_iav[i]);

 

 

 

T_iav[1]=funkf(0.0000001);

T_iav[N]=funkf(L-0.0000001);

 

}

 

 

void iavnii_metod_calc(){

time_iav=time_iav+tau_iav;

 

for(int i=1;i<=N+1;i++)

        TT_iav[i]=T_iav[i];

 

for(int i=2;i<=N-1;i++)

        T_iav[i]=TT_iav[i]+((lam*tau_iav)/(ro*c))*((TT_iav[i+1]-2*TT_iav[i]+TT_iav[i-1])/(h*h));

 

Form1->Series2->Clear();

for(int i=1;i<=N;i++)

        Form1->Series2->AddXY(h*i,T_iav[i]);

 

}

 

void __fastcall TForm1::FormCreate(TObject *Sender)

{

neiavnii_metod_start();

iavnii_metod_start();

}

 

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

void __fastcall TForm1::Button1Click(TObject *Sender)

{

time1=0;

neiavnii_metod_start();

iavnii_metod_start();

 

while (time1<t_max){

neiavnii_metod_calc();

iavnii_metod_calc();

}

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

}

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

 

 

void __fastcall TForm1::Button2Click(TObject *Sender)

{

time1=0;

neiavnii_metod_start();

iavnii_metod_start();

 

while (time1<TrackBar1->Position*tau){

neiavnii_metod_calc();

iavnii_metod_calc();

}

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

 

}

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

 

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

 

 

Вывод:

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

 

 

 

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

 

Вариант задания – метод Ньютона.

 

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

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

Пусть      − определена на отрезке     и дифференцируемая на нем вещественно-значная функция. Тогда формула итеративного исчисления приближений может быть выведена следующим образом:

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

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

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

Алгоритм:

1. Задается начальное приближение   .

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

Недостатки:

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

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

3. Если не существует вторая производная в точке корня, то скорость сходимости метода может быть заметно снижена.

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