CC BY
332
УДК 004.384:004.021
ОПТИМИЗАЦИЯ КОЭФФИЦИЕНТОВ ПЕРЕДАТОЧНОЙ ФУНКЦИИ, ПОЛУЧЕННОЙ МОДИФИЦИРОВАННЫМ МЕТОДОМ СИМОЮ ПО ЭКСПЕРИМЕНТАЛЬНО СНЯТОЙ ПЕРЕХОДНОЙ ХАРАКТЕРИСТИКЕ
А. Л. Рутковский, Л.И. Матвеева, Г.В. Козачек
Разработан программный комплекс и приведены результаты поиска оптимальных численных значений расчетных параметров передаточных функций, полученных модифицированным методом Симою, с целью минимизации ошибок аппроксимации исходных кривых разгона. Поставленная задача была решена методом наименьших квадратов с помощью функции MATLAB lsqnonlin (пакета Optimization Toolbox)
Ключевые слова: аппроксимация, передаточная функция, оптимизация
Анализ динамических характеристик
сложных технологических систем многомерных ТС составляет одну из ключевых задач проектирования алгоритмов управления и является весьма громоздкой вычислительной процедурой. При этом, ручной вариант расчета динамических характеристик возможен только в случае достаточно простых ТС.
Полное исследование динамических
характеристик ТС возможно только на математических моделях, адекватно описывающих временные связи параметров, законы движения и поведение моделируемого объекта в условиях действия возмущений.
Модифицированный метод Симою [1] включает последовательно выполняемые процедуры: 1)
выбор оптимальной структуры (формы) передаточной функции (из числа предлагаемых вариантов) исследуемого канала воздействия анализируемого динамического элемента и расчет ее параметров; 2) поиск оптимальных численных значений расчетных параметров полученных передаточных функций с целью минимизации ошибок аппроксимации исходных кривых разгона.
Алгоритмическую основу первой процедуры составляет один из распространенных методов расчета передаточных функций по кривым разгона (по переходным функциям) - метод Симою [2, 3], в основу которого положено предположение о возможности описания поведения исследуемого объекта (элемента) линейным дифференциальным уравнением с постоянными коэффициентами.
Для вывода основных соотношений, используемых при применении метода Симою, перейдем к операторному представлению динамики линейного объекта при воздействии единичного скачка на входе. Тогда изображение по Лапласу переходной характеристики (кривой разгона) представится в виде:
Рутковский Александр Леонидович - Северо-Кавказский ГМИ (ГТУ), д-р техн. наук, профессор, тел.(8672)57-43-08
Матвеева Людмила Ивановна - Северо-Кавказский ГМИ (ГТУ), канд. техн. наук, старший преподаватель, тел. 8906-494-33-75
Козачек Галина Владимировна - Северо-Кавказский ГМИ (ГТУ), соискатель, тел. 8-988-836-36-04 , e-mail -kzachek@rambler.ru
K
pF (p)
(1)
где W (p) =
bmpП + bm-1 p^ + .. + b1 p + 1 -
anpn + an-1 pn-1 +... + a1 p +1
передаточная функция линейного объекта. Для физически реализуемых систем п > т. Для общности рассмотрения принято, что п=т.
Введем в рассмотрение функцию
Хвъх ) = Хбых 0») - Хбых ) = К - Хвых (?) ,
изображение которой имеет следующий вид:
L [K - х„
х (f)] = ] [K-
X (f)]exp(- pf)df
K
-- У ы ( p)
F (p) - 1 K
F( p) p .
Разложим в ряд по степеням р изображение этой функции в окрестности точки р=0:
Ь [К - Хвых (0] = |[К - Хвых (0] ехр(-=|[К - Хвых (0] Л +
0 0
+ р\ [К - Хвых (0] ^ Ж + ... + р‘“^ [К - Хвъх (г)] & + ... = 2м, р ,
0 1 0 г- 1 =0
» 1 где М,. = Г [ К - Хвых (г)] — Лг - момент 1-го порядка, 1 1!
0
принимающий конечное значение на основании предположения того, что объект является звеном с самовыравниванием, т. е. отсутствуют
интегрирующие звенья.
Предположим, что кривая разгона снята с объекта с самовыравниванием, для которого справедливо
условие: Нш Х6ых (г) = К. Объекты без
самовыравнивания приводятся к объектам с самовыравниванием путем преобразования кривой разгона по известным методикам [2, 3].
Дальнейшие вычисления для них идентичны.
На основании полученных выражений,
аппроксимируем функцию ¥(р) в виде:
ад-^М— <2)
1 -2 КР"'
1=0 К
Определить коэффициенты
дифференциального уравнения при нулевых начальных условиях можно из условия выполнения равенства функции ¥(р) из (1) и правой части (2) для произвольных значений р, что приводит к известным выражениям метода моментов.
Однако наличие множителей (-г)1, входящих в подынтегральное выражение для вычисления моментов, уменьшает надежность определения коэффициентов передаточной функции в результате возрастания доли значений экспериментальных данных, измеренных при больших значениях г и имеющих наименьшую точность.
Для повышения надежности результатов метод моментов модифицируется следующим образом.
Предположим, что функция ¥(р) может быть представлена сходящимся рядом:
(3)
Тогда, приравнивая коэффициенты выражений (2) и (3) при равных степенях р, получим следующие рекуррентные соотношения:
51 = М 0; 5 2 = (1/К)51М0 + Мх;
8 3 = (1/К )[5 2М 0 + 51М1] + М 2;
Далее, используя формулу для моментов М1 и вводя новый временной масштаб в=г/Б1, можно исключить вычисление моментов и свести определение к однократному интегрированию (вычислению площадей) определенным образом модифицированной (приведенной к безразмерному виду [2, 3]) функции Хв^х(г):
»
81 = Г [К - (0№
0
» »
82 = 82 Г[К - ХвыхШа/К)-вв^[К - ХбЫхж2(еЛ0;
0 0
»
83 = 5? Г [К - (в)][(1/К) - 2в + (в2 / 2)]Лв =
0
»
= 82 Г [К - Хы (в)^_/ ;(в)йв;
0
84 = 84 Г [К - ХбЫх (в)][(83 / К813) - (82 / 82)в + в2 /2 - в3/6]Лв =
0 ■
= 84 Г [ К - Хвых (в)]Г4(в)йв
0
В общем случае произвольная 81 определяется следующим образом:
8, - 8'}[К - хв6„ (0)] 0
8,-1 + (-0)'-1 + (-0)і-2 + -'т3 5-1^(0 К8;-1 (і -1)! (і - 2)! £8;-1-у у!
^0
Приравнивая (1) и (3) и имея ввиду, что это равенство должно выполняться для любых р, получим:
81/ К = «1 - 61; 82 / К = а2 - 62 - (81 / К)Ь^
Для определения 2п коэффициентов а„ Ь1 совместно решается система уравнений,
включающая 2и-7
Ьп / й„ - х0 / К.
уравнений
(4)
С и™ :>....н
свод исходных данных
1. Шаг дискретизации (Ь)
2. Число интервалов дискретизации INI
3. Зю'ения (?) в N точках
4. Коэффициент усиления (К)
5. Ключ режима (Кл)
РАСЧЕТ времени чистого запаздывания (тг)
Смещение кривой БЛОК определения обве кга:
разгонана величину^ Х(і) = Х(і- Та ) - с самовыравниванием -без самоввіравнивания
Преобразование значений >0) ® безразтюрные величины
т
I
“ “С^^самиьыравнивание ”
ШОК подставления кривой разгона в виде разности двух кривых: 1. Эвена чистого запаздывания 2. Эвена с самовырвниванхем
^0)=1/№+1) ж
РАСЧЕТ
параметров
функции:
ЙТ» 0 ♦*)>)/ (гу+і>+і)
-----з-
РАСЧЕТ параметров функции: РГ(^)=(1+&1^ +ь?у
(г,р'+ту*г,Р*1)
Т
РАСЧЕТ параметров функции:
>) = (и-^р1+^1+&|;’)/
(1р+V* г^ + ^+О "Г
£
§____£
РАСЧЕТ РАСЧЕТ РАСЧЕТ
параметров функции: параметров функции: параметров функции:
1Г(р)=и ит» = н ИТ» = 1/
(г1р1 + ^+1) + 1) (і>- + іу+іУ+і>+ 1)
<р ________________
РАСЧЕТ параметров функции:
&м=а+ь,ру
(Т^+Т^+Т,р+1)
РАСЧЕТ параметров функции:
^С^)=(і +ъу+ь,ру {гу+тірі+ті?1+тір+1
£
РАСЧЕТ параметров функции:
^МІ+М' ,
(Г<р+7>+тіР^+71р+1)
Расчет параметров оригинала Х*(0 (обратное преобразование Лапласа)
ПЕЧАТЬ 1. Параметров W(p) 2.3начеши ошибки аппроксимации 3. Графики Х(1), Х*(1)
БЛОК оптимизации значении параметров передаточной
функции
Рис. 1. Блок-схема программного модуля расчета
параметров передаточной функции
и
,-1
т-1
Указанная процедура легко
алгоритмизируется.
На рис. 1 приведен алгоритм
модифицированного метода Симою.
Исходной информацией алгоритма является; N - число точек дискретизации интервала наблюдения: п - порядок линейного
дифференциального уравнения, описывающего
динамику объекта; к - шаг дискретизации
временной сетки; к - коэффициент передачи объекта; х1 - признак условия, что х6^1х(0)^0; {хбых(г), (1=1,... ,Щ}- значения абсцисс и ординат
кривой разгона (таблица значений кривой разгона). Результатом вычислений является совокупность коэффициентов дробно-рациональной
передаточной функции {ап, ап-1,.,а1, Ьп, Ьп-1,.,Ь1, к, Т, имеющей вид:
к(Ьпрп + Ьп ,рп-1 +... + Ь,р +1)
Ж(р) = -^-------п-^п~1------ехр(-т р).
апр + ап-1 р +... + а1 р +1
Программная реализация первой процедуры обеспечивает автоматизированный выбор оптимальной структуры (формы) передаточной функции (из числа предлагаемых вариантов) исследуемого канала воздействия анализируемого динамического элемента и расчет ее параметров с задаваемой точностью в условиях накладываемых ограничений по составу анализируемых видов (форм). Для наглядного представления исходной информации вводимые данные отображаются в графическом (построение исходной кривой разгона) и табличном (вводимые значения функции в различных точках) виде.
Результатом работы программы являются формы передаточных функций анализируемых каналов воздействия объекта, представленные в виде формул с численными значениями их параметров и графики, с изображением исходных и расчетных кривых разгона, а также численные значения оценок точности аппроксимации.
Программой предусмотрены два режима работы:
1) расчет параметров передаточной функции заданного вида (КЛ=1,...,7);
2) перебор возможных структур и расчет параметров передаточных функций (КЛ=0).
Максимальный суммарный порядок полиномов числителя и знаменателя в заложенных в программе базовых структурах (видах) передаточных функций равен 4. Программой предусмотрен предварительный анализ типа объекта (с самовыравниванием, без
самовыравнивания), преобразования кривых, расчет величины чистого запаздывания (тз) для приведения кривой разгона к стандартной форме, а также анализ наличия емкостного запаздывания (Т[), определяющего ветвь расчета. Для каждой анализируемой структуры передаточной функции проводится обратное преобразование Лапласа (разложение Хевисайда). Оценка точности полученных выражений передаточных функций программно осуществляется сравнением экспериментальных (полученных на динамических
объектах) с расчетными значениями динамической характеристики
1 N
Ay2 = N Жэ - V)2
п-\
где N - количество временных интервалов.
Практические расчеты передаточных функций по представленной методике показывают, что точность аппроксимации экспериментальной
кривой разгона не возрастает при вычислении площадей выше четвертого порядка. Однако,
остаточная квадратичная ошибка аппроксимации, при этом может быть достаточно большой. В связи с этим предложено использовать оценки
параметров передаточной функции как нулевое
приближение для решения задачи минимизации квадратичной ошибки аппроксимации.
Оптимизация передаточной функции,
определенной модифицированным методом Симою, заключается в минимизации суммы квадратов отклонений экспериментальных точек от решения дифференциального уравнения в тех же точках. Поставленная задача была решена методом наименьших квадратов.
Кривая разгона
[40 48 77 127 154 190 199 215 219 205 229 230 234 235 237 240
|«Г
О Data --------ОРТ
Рис. 2. Оптимизация коэффициентов передаточной функции
В качестве примера на рис. 2 показан результат аппроксимации данных
модифицированным методом Симою (кривая 8ішоіи) и результат оптимизации коэффициентов передаточной функции - кривая ОРТ.
Сформулируем задачу следующим образом: требуется определить значения
коэффициентов передаточной функции вида
W (p) =
коб
a3 p3 + a 2 p 2 + a1 p +1
ab a2, a3, такие, чтобы сумма квадратов отклонений экспериментальных точек от решения дифференциального уравнения в тех же точках была минимальна.
Функция MATLAB lsqnonlin (пакет
Optimization Toolbox) решает нелинейную задачу методом наименьших квадратов, включая задачу
нелинейного подбора данных с наименьшим числом итераций. x = lsqnonlin(fun,x0) начинает с точки х0 и находит минимум суммы квадратов описанной в fun функции. fun должна возвращать вектор значений, а не сумму квадратов значений. В данном алгоритме fun(x) неявно суммируется и возводится в квадрат.
До расчета значений f(x) (“суммы квадратов”), lsqnonlin требует задаваемой пользователем функции для расчета векторозначной функции.
fl(x)
F(x) =
f2(Х) f3(x)
Далее, используя задача оптимизации сформулирована как
векторную терминологию, может быть заново
тіп-іЦ F (X)
Поскольку lsqnonlin предполагает, что сумма квадратов явно не формируется в пользовательской функции, то передаваемая в lsqnonlin функция должна взамен этого вычислять векторозначную функцию, где х есть вектор и F(x) есть функция, которая возвращает значение вектора.
Вектор начальных значений а0 содержит подлежащие оптимизации значения а1, а2, а3, функция fun реализована как разность
экспериментальных данных и решения
дифференциального уравнения. function f=myf(a,Y) global s1 kob t a0 n dx j=1:n;
W=tf([kob*a(3),kob],[a(2),a(1),1]);
x=dx*ones(1,n);
Y=s1(1)+lsim(W,x,t);
f=s1(j)-Y;
W - передаточная функция;
У - решение дифференциального уравнения в фиксированных точках;
f - разность экспериментальных данных и решения дифференциального уравнения;
81 - массив экспериментальных данных.
В приведенном примере исходные значения параметров передаточной функции а1=8.6089, а2=18.6371, а3=3.2802 были определены
модифицированным методом Симою и являлись начальными значениями для lsqnonlin, которая возвратила вектор значений а1 =8.7759, а2= 16.3827, а3= 28.9226.
Аналогично производится оптимизация для любой другой структуры передаточной функции.
Литература
1. G.G. Arunyants, A.L. Rutkovskii, Z.G. Salikhov, and D.N. Stolbovskii, Computation of Dynamic Characteristics of Control Systems An Effectiveness Enhancement Method, Automation and Control, Vol. 66, № 4, 2005, pp. 562-569
2. Симою М.П., Определение коэффициентов
передаточных функций линеаризованных звеньев и систем авторегулирования, «Автоматика и
телемеханика», т. XVIII, 1957, № 6
3. Стефани Е.П. Основы расчета настройки регуляторов теплоэнергетических процессов. М-Л: Государственное энергетическое издательство, 1960. С. 328.
4. А.Г.Трифонов. А.Г."Постановка задачи оптимизации и численные методы ее решения" сайт exponenta.ru pаздел "Математика\Optimization Toolbox
2
2
Северо-Кавказский горно-металлургический институт (Государственный технологический университет)
OPTIMIZATION OF FACTORS OF TRANSMISSION FUNCTION,
GOT BY MODIFIED METHOD BY SIMOYU ACCORDING TO THE EXPERIMENTAL
SKIM CONNECTION FEATURE
A.L. Rutkovskiy, L.I. Matveyeva, G.V. Kozachek
The program complex has been designed and the results of searching of optimum numerical importances accounting parameter transmission functions, got by modified method Simoyu for the reason of minimization mistake to approximations of the source curve of runway has been brought. The put problem is solved by method the least square by means of functions MATLAB
Key words: approximation, transmission function, optimization
