Известная функция 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)):

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

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

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

Страницы: 1, 2, 3