Построение математических моделей методом идентификации (стр. 4 из 5)
Продолжим этот процесс до 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. Численные процедуры оценивания параметров нелинейных регрессионных моделей
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. Разработка аналитических моделей объектов автоматизации. Линеаризация моделей
5.1 Постановка задачи
Разработайте математическую модель заданной динамической системы. Построить аналитическую модель емкости с одним притоком и двумя стоками (например двухручьевая МНЛЗ) и провести ее линеаризацию. Математическую модель заданного объекта представить в виде программы на ЭВМ, рассчитать кривые исходной и линеаризованной модели отклика модели на ступенчатое и импульсное воздействие. Предложить методику идентификации параметров модели. Какие эксперименты будет необходимо провести на реальном объекте? Оцените необходимый объем и форму представления результатов.
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 – уровни металла в сталеразливочном и промежуточном ковшах;