авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 |   ...   | 4 | 5 || 7 | 8 |   ...   | 9 |

«Е.А. Новиков, Ю.В. Шорников КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ЖЕСТКИХ ГИБРИДНЫХ СИСТЕМ Е.А. Новиков, Ю.В. Шорников КОМПЬЮТЕРНОЕ ...»

-- [ Страница 6 ] --

Пусть отображение C ограничено классом систем обыкновенных дифференциальных уравнений (ОДУ) в форме Коши с запаздывающим аргументом, т. е.

x = f [ x(t ), x(t ), t ], x(t0 ) = x0, (9.1) где x(t ) R n – вектор состояния;

x(t ) = (t ) при t [t0, t0 ) ;

t – не зависимая переменная;

(t ) = ( 1(t ),..., s (t ) ), s n – s-мерная вектор Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ функция запаздывания;

= {1,..., s }T – вектор чистых запаздываний;

f : R R n R n – нелинейная вектор-функция, удовлетворяющая ус ловию Липшица;

x0 R n – вектор начальных условий.

Дискретными называются такие объекты, состояние которых опре деляется вектором d W r, W = R B, – множество неотрица тельных целых;

B = {false, true} – множество булевых величин.

Поведение дискретной системы задается отображением D : R W r W r, которое определяет причинно-следственную цепоч ку событий:

E0 = (t0, x0 ) E1 = (t1, x1 )... E j = (t j, x j )....

В общем случае поведение С гибридной системы зависит от сово купности поведений {c j, j = 1,..., m}. Качественные изменения поведе ний непрерывной системы происходят мгновенно в некоторые момен t = t* ты времени для всей совокупности поведений. Область j существования поведений {c j, j = 1,..., m} в расширенном фазовом (t, x) R n+ пространстве задается логическим предикатом pr j (t, x) B, который принимает значение false внутри области суще ствования c j и значение true – вне ее. Таким образом, момент време ни t = t * наступления очередного события E j определяется моментом j изменения pr j (t, x) B от значения false на значение true. Несмотря на то что в любой физической системе время t R априори является величиной непрерывной, можно говорить о дискретных моментах вре мени t * как подмножестве значений непрерывного времени. Такие j моменты принято называть временной щелью (time gap), и в этом слу чае рассматривается так называемое гибридное время (t, k ) Th : R.

Гибридное время имеет вполне определенный физический смысл. Рас смотрим в связи с этим эффект Зенона.

9.2. Эффект Зенона 9.2. ЭФФЕКТ ЗЕНОНА Пусть имеются два бака одинакового размера (рис. 9.1), из которых вытекает жидкость со скоростью v1 и v2 соответственно. Для каждого бака существует свой предельный уровень жидкости r1 и r2, после чего следует доливать жидкость в баки. Кран для доливки жидкости один, из него жидкость поступает в любой бак со скоростью w. Кон троллер в зависимости от предельных уровней r1 и r2 мгновенно пе реключает поток w.

w w 1 x r Рис. 9.1. К эффекту Зенона Предположим, что жидкость начинает доливаться в тот момент, ко гда она в первом баке достигла критического уровня, а во втором его превышает. Доливка жидкости в первый бак будет продолжаться до тех пор, пока уровень жидкости во втором баке не достигнет критиче ского уровня, после чего начинает пополняться второй бак. Гибридная система будет иметь два состояния c1 и c2 с уравнениями dx1 dx c1: x2 r2, = w v1, = v2 ;

dt dt dx1 dx c2 : x1 r1, = v1, = w v2.

dt dt Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ Положим r1 = 0 и r2 = 0. Начальным будем считать состояние, ко гда x1 = 0 и x2 = h 0, т. е. в момент времени 0 = 0 при x1 ( 0 ) = 0 и x2 ( 0 ) = h начинается заполнение первого бака. В момент времени 1 = x2 (0 ) v2 = h v x2 (1 ) = 0. За время уровень во втором баке станет равным 1 = 1 0 первый бак заполнится до уровня x1 ( 1 ) = ( w v1 )1 = ( w v1 ) h v2.

Произойдет переключение, и начнет заполняться второй бак. В момент h w v x ( ) h h 2 = 1 + 1 1 = + ( w v1 ) = 1 + v1 v2 v2 v1 v2 v второй бак за время ( w v1 )h 2 = 2 1 = v1v будет заполнен до уровня ( w v2 )( w v1 )h x2 (2 ) = ( w v2 )2 =.

v1v Первый бак станет пустым, и снова произойдет переключение.

В момент следующего переключения x ( ) 3 = 2 + 2 2 = v h w v1 h = 1 + + ( w v1 )( w v2 ) = v2 v1 v1v2 v h w v1 w v = 1 + 1 + v2 v1 v в первом баке будет уровень ( ) x1 ( 3 ) = ( w v1 )3 = ( w v1 ) 2 ( w v2 )h v1v2.

9.2. Эффект Зенона Введем обозначения:

1 = ( w v1 ) v1 и 2 = ( w v2 ) v2.

Очевидно, что имеет место h h (1 + 1 ) j 1 (1 + 2 ) j 1, 2 j = (1 + 1 ) j (1 + 2 ) j 1, 2 j 1 = v2 v где 2 j 1 = h1j 1 21 v2, 2 j = h1j 2 1 v2, j = 1, 2,….

j j Суммируя по j, имеем h h 1j 21 + 1j 121 = j j j = 2 j + 2 j 1 = v v j =1 j =1 j =1 j =1 2 j =1 h (1 + 1 ) ( 12 ).

= v2 j = Допустим, что выполняется неравенство max {v1, v2 } w v1 + v2.

Тогда 12 1, и в итоге получим v1v (12 ) j = 1 =.

w(v1 + v2 w) j =0 Окончательно имеем h( w v1 ) hv1 h k = w(v + =.

1 + v2 w) w(v1 + v2 w) v1 + v2 w k = Таким образом, объединение всех интервалов последовательности, образующей гибридное время, не совпадает с положительной вещест венной осью, т. е.

k = const.

k = Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ Последнее условие можно интерпретировать как переключение с бесконечной скоростью в окрестности некоторой временной точки.

Такую ситуацию будем называть эффектом Зенона.

Введение гибридного времени Th позволяет рассматривать гло бальное поведение гибридной системы как совокупность локальных поведений {c j, j = 1, m}, причем решения в моменты переключения () x j t* являются начальными условиями x 0+1 для нового поведения j j с j +1. Глобальное поведение гибридной системы – как бы «склейка»

локальных поведений. Заметим, что в теории дифференциальных уравнений говорят о «припасовывании» решений вместо термина «склейка» локальных поведений в гибридных системах. Отметим так ( ) ( ) x j t* 0 x j t* + же, что в общем случае имеет место j j j J 1... m, поскольку не всегда c j = c j +1.

9.3. РЕЖИМЫ И СОБЫТИЯ Определение 9.1. Режимом ГС будем называть кортеж nj nj x0j j j j, j {1... nm }, f : R R R j, x, f,, где режимное пове nj характеризуется начальным значением x j (t0 ) = x0j, дение x j (t ) R а режимная функция f j (t, x j ) удовлетворяет условию Липшица для n всех (t, x j ) R R j.

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

Определение 9.2. Непрерывно дифференцируемая функция g (t, x) : R R n R s, s = 1, 2,..., называется событийной функцией и характеризуется предикатом pr : g (t, x) 0, pr B = {false, true}.

Событийная режимная функция g j (t, x) 0 ведет себя таким об разом, что соответствующий предикат pr j : g j (t, x) 0 режимного 9.3. Режимы и события поведения является истинным ( pr j = true) на всем полуинтервале ре ) жимного решения t 0, t * [t0, tk ].

j j Пример 9.1. Рассматривается элементарная классическая гибрид ная система – абсолютно упругий отскок мячика от абсолютно твердой горизонтальной поверхности. Начальное положение мячика над по верхностью отскока y = H. Описание системы с учетом только верти кальной составляющей скорости v y и расстояния y до поверхности отскока описывается уравнениями y = v y, v = q, y y (0) = H, v y (0) = 0, if ( y 0) and (v y 0) then v y = v y, где q – ускорение свободного падения. Это типичный случай гибрид ной системы с мгновенным изменением значений фазовых координат по условию достижения мячиком поверхности отскока, когда y 0 и vy 0.

Режим ГС может быть ограничен одним или несколькими усло виями. В данном примере режим полета до поверхности отскока соот ветствует одновременному выполнению двух условий:

1) вертикальная составляющей скорости положительна, т. е.

g1 (v y, t ) = v y 0 ;

2) расстояние до поверхности отскока – также величина положи тельная, т. е. g 2 ( y, t ) = y 0.

Событийная функция g ( t, x ) в общем случае нелинейна. Введени ем дополнительной переменной z = g (t, x) получим в качестве ограни чений новую линейную событийную функцию x = f ( x) x = f ( x) g g z = f ( x) +. (9.2) g (t, x) 0 x t z 0 Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ Можно перейти от z к расширенному фазовому вектору x по правилу zi = xn +i, 1 i n. Тогда исходная система перепишется в виде задачи с линейной событийной функцией, т. е.

x = f ( x), xi 0, i = n + 1,..., 2n, где x R 2n, f : R R 2n R 2n, причем f n +i = gi (t, x) t, 1 i n.

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

Определение 9.3. Границей режима гибридной системы называ ется некоторая граница области G R n, на которой для событий ной функции имеет место g (t, x) = 0, при этом t = t *.

Отметим, что точность поиска точки t * переключения режимов в гибридных системах существенно зависит от вида событийной функ ции g (t, x), линеаризация которой не только увеличивает размерность задачи, как это следует из (9.2), но и приводит к понижению точности определения точки переключения режимов, если не использовать спе циальный алгоритм.

Определение 9.4. Событием гибридной системы будем называть такое состояние ГС в пространстве и времени (t, x) R n+1, когда g (t, x) : R R n R s достигает границы режима.

События делят на три категории: односторонние, двусторонние и критичные к точности события. Так как фазовые траектории вычис ляются с помощью арифметики с конечной точностью, то такое деле ние необходимо, потому что невозможно точно определить моменты переключения. Определение момента времени t = t *, когда событийная функция g ( t, x(t ) ) = 0, является достаточно сложной задачей коррект ного обнаружения событий. Ситуации, когда существуют особые точ ки или когда физический смысл проблемы диктует такое условие, что фазовая траектория никогда не должна пересекать поверхность (грани цу) события, относятся к категории односторонних событий.

Определение 9.5. Односторонними событиями гибридных систем будем называть события, для которых имеет место x 0.

9.4. Локальное и глобальное поведение Иллюстрацией систем с односторонними событиями являются, на пример, модели реакций химической кинетики. Динамика концентра ции компонентов с учетом их физического смысла может быть только величиной неотрицательной. Идентификация процесса химической кинетики выполняется только с ограничениями g ( t, x(t ) ) 0, x 0, учитывающими односторонние события. Пример прыгающего мячика также является типичным примером гибридной системы с односто ронними событиями. Именно такие события представляют наиболь ший практический интерес и поэтому будут исследоваться в дальней шем. Гибридные системы в данном случае относят к дискретно непрерывным, причем их непрерывное поведение в основном анализе рассматривается как численное решение задачи Коши. В предыдущих главах предложены оригинальные численные схемы для решения же стких систем дифференциальных уравнений. Актуальность разработки надежных эффективных численных методов решения жестких задач распространяется и на гибридные системы.

Определение 9.6. Жесткой гибридной системой будем называть такую систему, в которой хотя бы один из режимов является жестким.

Жесткость режима в данном случае рассматривается в смысле оп ределения 2.1. Из введенных формальных определений следует, что динамические системы являются подклассом гибридных систем или режимами ГС. В связи с этим локальное поведение гибридных систем можно рассматривать как поведение динамических систем в расши ренном пространстве состояний.

9.4. ЛОКАЛЬНОЕ И ГЛОБАЛЬНОЕ ПОВЕДЕНИЕ Локальное поведение гибридной системы характеризуется единст венным режимным поведением, полученным при решении задачи Ко ши с ограничениями:

y = f ( t, y ), y (t0 ) = y0, t [t0, tk ], (9.3) pr : g (t, y ) 0, где y R N – вектор состояния;

f : R R N R N – нелинейная век тор-функция, удовлетворяющая условию Липшица;

y0 R N – вектор начальных условий.

Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ Здесь и ниже, в отличие от (9.1), будем рассматривать функцию f без отклонения аргумента, поскольку (9.1) легко сводится к (9.3) ап проксимацией Падэ с повышением порядка системы до N, N n, где n – исходный порядок системы уравнений (9.1). Процедура замены (9.1) системой (9.3) состоит в следующем. В структурных схемах, со ответствующих (9.1), звено запаздывания можно представить переда точной функцией в виде отношения двух полиномов следующего вида:

L M H ( e p ) = Bi pi Aj p j, i =0 j = где ( L + M i )! (1)i i M d, i = 0, L, L M ;

p= ;

Bi = ( L i )!

dt i!

( L + M j )! M ! j M Aj =, j = 0, M.

( M j )! j !

L!

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

y u Рис. 9.2. Структурная схема, соответствующая канонической форме аппроксимации Падэ Таким образом, структурная схема, замещающая рассматриваемую передаточную функцию запаздывания, состоит из nk интеграторов.

9.4. Локальное и глобальное поведение В обратной связи находятся коэффициенты знаменателя исходной пе редаточной функции, а в прямой связи – коэффициенты полинома чис лителя передаточной функции. В результате замещения порядок сис темы увеличивается до N = n + nk, где n – порядок исходной системы уравнений (9.1), nk – количество звеньев запаздывания в (9.1). Из оп ределений 9.1 и 9.2 следует, что поскольку ограничения в (9.3) всегда выполняются, их можно опускать, не нарушая общности. Тогда (9.3) является классической задачей Коши (2.1), на которую в дальнейшем и будем ссылаться, если речь идет о режимном или локальном поведе нии гибридной системы.

Режимные поведения некоторых гибридных систем описываются дифференциально-алгебраическими уравнениями (ДАУ). Примеры таких систем рассмотрены ниже. Расширим класс задач (9.3) введени ем алгебраических уравнений, при этом ограничимся алгебраическими уравнениями, разрешенными относительно алгебраических перемен ных z. Тогда (9.3) примет вид y = f (t, z, y ), z = (t, z, y ), y (t0 ) = y0, z (t0 ) = z0, (9.4) pr : g (t, y ) 0, t [t0, tk ], где y, z R N, f : R R N z R N R N, : R R N R N z R N z.

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

Семантика ограничений g ( x, t ) 0 связана с событийными функ циями гибридных систем. Временной интервал существования локаль ного поведения ГС определятся временным отрезком j t0j, tkj [t0, tk ], j = 1, m, на котором условие g ( y, t ) 0 истинно. Говорят, что система на этом отрезке находится в некотором локальном состоя нии. Мгновенный переход ГС в новое локальное состояние с новым Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ поведением, определенном решением (9.3), происходит при наруше нии условия ограничения режима ГС, т. е. когда предикат принимает значение pr = false. Считается, что смена локальных состояний проис ходит мгновенно в некоторый момент времени t * [t0, tk ].

Глобальное поведение ГС характеризуется совокупностью согласо ванных режимных поведений, полученных на множестве решений за дачи Коши с ограничениями вида y j = f j (t, y j ), y j (t0 ) = y0, (9.5) j t [t0, tk ], pr j : g j (t, y ) 0, где Nj Nj Nj Nj sj y j R, f j :RR R, g j :R R R, s j = 1, 2,…, 1 j m.

В настоящее время актуальны проблемы поиска решения событий но-непрерывных систем вида (9.5). Попытки применить классические подходы к анализу гибридных систем пока дают весьма ограниченные результаты. В связи с этим исследования таких систем необходимо проводить методом компьютерного анализа с использованием инстру ментальных средств.

Итак, в общем случае можно сказать, что под гибридной понимает ся система, поведение которой характеризуется совокупностью непре рывных локальных поведений при мгновенных дискретных переходах из одного локального состояния в другое. Совокупность дискретных переходов обусловлена непустым множеством событийных функций g = {g j }, 1 j m.

Время существования режима ГС определяется интервалом между двумя соседними событиями. Условия Липшица нарушаются в момен ты времени t * наступления очередного события c j. В точке разрыва j не существует производной f (t, y ), что вызывает нарушение условия Липшица.

9.5. Разрывы 9.5. РАЗРЫВЫ Дифференциальные уравнения с разрывной правой частью рас сматривались и ранее в связи с тем, что такие системы часто встреча ются в механике, электротехнике, теории оптимального управления и других областях науки. Источником задач с разрывной правой частью в механике служат, например, механические системы с сухим трением.

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

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

Нелинейные системы автоматического управления с переключе ниями, а также гибридные системы с мгновенными изменениями ло кальных состояний в результате изменений начальных условий и пра вой части описываются в классе уравнений (9.3) с разрывными по x функциями f (t, x). Остановимся подробнее на нарушении условий Липшица. Рассмотрим систему y = f (t, y ), y (t0 ) = y0 (9.6) и начнем с более простого случая, когда правая часть дифференциаль ного уравнения в (9.6) непрерывно зависит от y и может быть разрыв ной только по t. В случае непрерывной функции f задача Коши для уравнения (9.6) эквивалентна интегральному уравнению t f ( s, x( s) ) ds.

y (t ) = y0 + (9.7) t Один из наиболее известных и распространенных наборов условий на правую часть уравнения (9.6), при которых такой подход к расши рению понятия решения оказывается содержательным в математиче Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ ском плане и полезным и адекватным в приложениях, выглядит сле дующим образом. Предполагается, что выполнятся следующие усло вия Каратеодори:

1) функция x f (t, y ) непрерывна почти при всех t ;

2) функция t f (t, y ) измерима при каждом x ;

3) существует локально суммируемая функция u : R R такая, что при каждом фиксированном y почти при всех t выполняется не равенство f (t, y ) u (t ), где || || – некоторая норма в R N.

Эти условия, в частности, гарантируют, что если функция t y (t ) измерима, то функция t f [t, y (t )] и решением Каратеодори задачи (9.5) на промежутке [t0 T, t0 + T ] называют абсолютно непрерывную функцию, удовлетворяющую начальным условиям задачи и обращаю щую почти всюду на [t0 T, t0 + T ] уравнение (9.7) в тождество.

Распространение теоремы существования решения на дифференци альные уравнения с разрывными правыми частями связано с обобще нием понятия решения, поскольку f (t, y ) не всюду может иметь про изводную в связи с разрывами. Согласно теореме Каратеодори, решением задачи (9.6) является решение интегрального уравнения (9.7) в смысле Лебега.

Нелинейные системы автоматического управления с переключе ниями и гибридные системы описываются в классе уравнений (9.6) с разрывными по x функциями f (t, y ). Достаточное условие единствен ности решения (9.6) при t t0 заключается в выполнении неравенства ( f (t, y1 ) f (t, y2 ) )( y1 y2 ) k (t ) y1 y2, где k (t ) – интегрируемая в области D R N +1 вектор-функция, (t, y1 ), (t, y2 ) D R N +1.

В этом случае для любого начального условия y (t0 ) = y0, (t0, y0 ) D при t t0 может существовать не более одного решения задачи (9.5).

9.6. Анализ событийно-непрерывных систем 9.6. АНАЛИЗ СОБЫТИЙНО-НЕПРЕРЫВНЫХ СИСТЕМ Приведем примеры событийно-непрерывных систем. В первом примере ставится вопрос о существовании автоколебаний в нелиней ных моделях второго порядка, в частности в модели часов с различны ми типами мгновенного удара спускового механизма:

dx + f0, 0, dt dx +x= 2 dx dt f, 0, dt где f 0 – постоянная сила трения, отнесенная к единице массы. Удар подчиняется либо предположению о постоянстве количества движения v1 v0 = const, либо предположению, что кинетическая энергия систе мы при ударе изменяется на одну и ту же величину:

2 mv1 2 mv0 2 = const.

Здесь v1 и v0 – скорости после и до удара. Введением новой фазовой переменной y = dx dt модель легко приводится к форме (9.5), т. е.

c1 : x = y, y = f 0 x;

pr : g1 = y 0;

c2 : x = y, y = f 0 x;

pr : g 2 = y 0.

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

Следующий пример – модель автопилота с релейным сервомото ром описывается уравнениями d 2 d 2 = M h () + q, dt dt d = +, () = sign(), dt Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ где q, M, h и – константы. Эта система с так называемым релей ным управлением является гибридной системой с двумя локальными поведениями: c1, если g () = 0, и c2, если g ( ) = 0. Вариаци ей параметра можно улучшать управление автопилотом. Как видно из рис. 9.3, при увеличении переходной процесс заканчивается бы стрее, чем без управления по производной, когда = 0. Однако с дальнейшим увеличением возникает эффект Зенона (рис. 9.3, б).

Этого эффекта можно избежать заменой разрывной функции () = sign() гладкой функцией () = 2arctg(K ).

0, 0, –0, –0, –0, –0, –0,0352 0 0,0352 0,0704 0,106 –0,0315 0,063 0, 0 0, а б Рис. 9.3. Улучшение управления автопилотом за счет увеличения значения коэффициента : 1 Тогда поведение ГС вырождается в поведение обычной динамической системы = x, x = Mx 2arctg(K ), = + x с результатами, приведенными на рис. 9.4.

Исследование показывает, что при использовании функции () = sign() появляется эффект Зенона, а для функции () = 2arctg ( K ) начинается движение к точке стабилизации, как это видно из рис. 9.4. Данное явление называется явлением скольже ния. При использовании модели с разрывной функцией () = sign() следовало бы обнаружить эту прямую и доопределить уравнения 9.6. Анализ событийно-непрерывных систем движения. При использовании гладкой функции () = 2arctg ( K ) это произошло автоматически.

х 0, –0, –0, –0, 0, –0,0319 0, 0 0, Рис. 9.4. Аппроксимация функции () = 2 arctg ( K ) Явление скольжения рассматривается и в работах Е.С. Емельянова при решении задач стабилизации. Аналитические исследования таких систем иллюстрируются на примере линейной системы второго порядка:

x1 = x2, x2 = a2 x2 a1 x1 bu, где b 0, c 0, – константы,, x1 ( x2 cx1 ) 0, u = x1, =, x1 ( x2 cx1 ) 0.

В частности, показано, что введением управления = система ока зывается колебательно неустойчивой, а при = – апериодически неустойчивой.

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

Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ x1 = x2, x2 = a2 x2 a1 x1 bx1u, u=, где || || – некоторая норма в R n. В качестве нормы функции управле ния предлагается считать следующую норму:

+ ( ) = arctg ksigma x1 ( x2 + cx1 ) +.

16,4 х 14, х2 8, 7, –8, 14,4 х х –16, 27, 33,8 –13,8 0 13, –8,45 0 8,45 16,9 25,4 –27, а б Рис. 9.5. Явления скольжения (а) и постоянного переключения (б) Полученные при этом результаты анализа с параметрами = 1 и = 1 приведены на рис. 9.5. Вариацией параметра c система с пере менной структурой демонстрирует себя как систему с постоянным пе реключением, иногда в системе наблюдается устойчивое движение с явлением скольжения.

9.7. КЛАССИФИКАЦИЯ СОБЫТИЙ При возникновении событий в системах, где длительное поведение описано с помощью обыкновенных дифференциальных уравнений, происходит мгновенное изменение локальных состояний или режимов гибридных систем. Переключение режимов гибридных систем приво дит к разрывам в задаче (9.5). Данные разрывы обусловлены скачкооб разным изменением начальных условий, параметров правой части, из менением правой части без изменения и с изменением состава фазовых переменных и др. Рассмотрим решение перечисленных ситуаций в свя зи с наступлением событий в ГС.

9.7. Классификация событий Изменение начальных условий Скачкообразное изменение значений некоторых переменных {xi | i J 1, n} при x(t * 0) x(t * + 0) можно рассматривать как ре зультат действия некоторой последовательности операторов присваи вания xi := xi* | i J 1, n, выполняемых мгновенно в момент переключения t *. В этот момент происходит припасовывание начальных условий. Иллюстрацией к этому типу событий является рассмотренный выше пример с абсолют но упругим отскоком падающего тела. При касании тела поверхности отскока, когда расстояние до поверхности отскока y 0 и вертикаль ная составляющая скорости v y 0, предикат ограничений pr : g ( y, v y ) = { g1 ( y ) = y 0;

g 2 (v y ) = v y 0} принимает значение pr = false, а вертикальной составляющей скоро сти разрешается мгновенно сменить знак на противоположный, т. е.

v y = v y. После этого снова начинается интегрирование исходной системы с новыми начальными условиями. Таким образом, поведение гибридной системы может быть вполне удовлетворительно описано последовательностью из припасованных решений дифференциальных уравнений на отдельных участках решения. На стыках длительных по летов можно пренебречь длительностью реального отскока, превратив временной промежуток в точку – достаточно считать, что скорость мгновенно меняет знак на противоположный.

Изменение параметров правой части Представим (9.5) с изменяемыми параметрами R N в следую щем виде:

y = f ( y, t, ), y (t0 ) = x0, t [t0, tk ], = j, 1 j N, Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ и рассмотрим изменение значений параметров R N двумя спосо бами:

1) получение новой задачи введением новых фазовых переменных;

2) введение бинарных компонент в правую часть исходной задачи (9.5).

Расширение вектора фазовых переменных Введем набор параметров R N в состав фазовых переменных.

Тогда, если положить = const, новую модель гибридной системы можно представить в виде y = f ( y, t, ), y (t0 ) = y0, = 0, (t0 ) = 0, N j = R n, t [t0, tk ] j, j = с расширенным вектором состояний ( y, ) R N + N. Затем решение ищется припасовыванием решений классической задачи Коши на уча ) стках t0j, t * с неизменными параметрами и мгновенно меняющимися j начальными условиями в моменты времени t *. Для иллюстрации ис j пользования приема расширения вектора R фазовых переменных рассмотрим пример гибридной системы, в которой событием является изменение параметров правой R С части.

Е Пример 9.2. Пусть в электрической це пи (рис. 9.6) в зависимости от положения ключа меняется значение сопротивления L = R.

Рис. 9.6. Элементарная элек Поведение этой цепи описывается трическая цепь с ключом уравнением t di + iR + id = E, L dt c 9.7. Классификация событий и если R меняется периодически R (t ) = R(t + T ), например, при T = по формуле R1, t = 0,3,…;

R= R2, t = 1, 2,…, то получим систему с мгновенным изменением параметров в моменты времени t * = 1, 2,.... Преобразуем ее к виду t di du + iR + u = E, u = id, = i, L dt dt c c и дополним новыми дифференциальными уравнениями di 1 du 1 dR = ( (iR + u ) + E ), = i, = 0.

dt L dt c dt Последняя система будет эквивалентна предыдущей, если рассматри вать ее как набор систем, у которых периодически меняются началь ные условия последнего уравнения в точках t * = 1, 2,.... В эти моменты времени мгновенного изменения начальных условий фазовые пере менные i и u «склеиваются», образуя непрерывные функции, а R ме няется скачком, т. е. ведет себя как кусочно-непрерывная функция.

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

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

g (t ) = sign r, где || || – некоторая норма в линейном вещественном пространстве R, а функция sign(r ) задается формулой Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ 1, r 0;

sign(r ) = 0, r = 0.

Тогда значения моментов t * наступления событий переключения эк вивалентны нулям или единицам событийной функции g ( t ). Очевид но, что событийная функция t t g (t ) = sign T T удовлетворяет условиям сформулированной задачи периодических пе реключений параметров правой части ОДУ в дискретные моменты времени t * с периодом T, где [] – целая часть числа. Тогда обобщен ная система дифференциальных уравнений рассмотренной задачи бу дет иметь вид di 1 du = ( (iR + u ) + E ), = i, dt L dt c R = R1 g (t ) + R2 R2 g (t ), причем параметры правой части будут меняться автоматически по за кону изменения событийной функции. Событийная функция в этом случае ведет себя как бинарная компонента, принимающая значение нуля или единицы, и она же управляет параметрами правой части сис темы ОДУ.

Для рассмотренного примера с вычислительной точки зрения про ще использовать прием с автоматическим изменением параметров вве дением бинарных компонент в правую часть. Это можно сделать, если удается подобрать простой закон изменения параметров правой части, не включая в программную модель условные операторы для изменения значений параметров. Кроме того, использование приема введения до полнительных фазовых переменных увеличивает размерность исход ной задачи. Прием введения бинарных компонент предпочтительнее еще и потому, что нет необходимости при реализации пользоваться условным оператором для переопределения начальных условий, как в примере 9.1. Переопределение начальных условий или новых парамет ров правой части будет происходить автоматически в моменты t * че рез событийную функцию.

9.7. Классификация событий Изменение вида правой части без изменения состава фазовых переменных В этом случае происходит изменение отображения C. Примером такой гибридной системы является падение мячика на вертикально стоящую пружину, после которого в уравнении движения мячика по является сила реакции пружины. Пусть мячик падает на свободный конец невесомой пружины с жесткостью K и длиной H s H, закреп ленной вертикально на плоскости. Движение мячика при y H s зада ется системой уравнений свободного падения (см. пример 9.1) и при g ( y ) = y H s – системой уравнений dv y dy = K ( H s y ) q.

= vy, dt dt Результаты численного моделирования представлены на рис. 9.7.

v(y) 9, 6, 3, –3, –6, y –9, 8, 7,33 10, 5, 1,47 2,93 4, Рис. 9.7. Мячик, падающий на пружину. Фазовая диаграмма Изменение правой части режима ГС с изменением состава фазовых переменных Рассмотрим последний случай, когда наступление событий в гиб ридной системе приводит к качественным изменениям локального по ведения.

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

Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ Состояние маятника в режиме колебаний определяется уравнениями y =, = g sin() L, x = L sin(), y = L cos(), x v x = cos(), v y = sin(), (9.8) g () = max 0, g (0) = 0, (0) = 0.

Пусть в некоторый момент t* крепле Рис. 9.8. Отрывающийся ние шарика к стержню разрушается и да маятник лее шарик продолжает свое независимое от стержня движение. Динамика маятника до отрыва определяется двумя дифференциальными уравнениями и ограничением на событий ную функцию, определяющую новое событие – момент отрыва. Дви жение шарика после отрыва задается новой системой уравнений с но вым составом фазовых переменных и, как в примере 9.1, с учетом горизонтальной составляющей может быть представлено в виде dv y dx dy = q.

= vx, =vy, (9.9) dt dt dt 0 y(x) –0, –0, –0, –0, –0, –0,95 x 0,402 0,803 1,2 1, –0,803 –0, Рис. 9.9. Фазовая диаграмма отрывающегося маятника Начальные значения для новой системы уравнений, описывающей свободный полет, вычисляются в момент t* по формулам 9.8. Инструментально-ориентированный анализ x(t * ) = L sin(), y (t * ) = L cos(), v x (t * ) = cos(), v y (t * ) = sin().

На рис. 9.9 показана траектория движения маятника для значений па раметров 0 = 2, 0 = 0 и max = 4.

9.8. ИНСТРУМЕНТАЛЬНО-ОРИЕНТИРОВАННЫЙ АНАЛИЗ Один и тот же класс задач можно изучать с помощью различных подходов: с помощью аналитического способа для анализа систем автоматического управления и численного подхода для анализа гиб ридных систем. В настоящее время актуальна проблема поиска общего решения событийно-непрерывных систем вида (9.5).

При численном анализе гибридных систем целесообразнее исполь зовать явные одношаговые схемы по следующим причинам. В гибрид ных системах, как и в обычных динамических системах, для поиска решения xs неявными методами приходится прибегать к итерацион ным процедурам, например методу Ньютона–Рафсона с дорогостоя щей вычислительной процедурой вычисления и обращения матрицы Якоби. Для k -шаговых схем типа Адамса необходимо знать и помнить значения фазовых координат xs i и их производные xs i на k преды дущих шагах. Однако в момент запуска и после точек разрыва при мгновенном переходе в новые локальные состояния при определении этих величин возникают некоторые трудности. Кроме того, использо вание явных одношаговых численных алгоритмов менее опасно в ус ловиях односторонних событий гибридной системы. Следует отме тить, что в отличие от L -устойчивых методов явные схемы имеют ограниченные области устойчивости, что мешает увеличению шага интегрирования на участках установления, где шаг ограничен критери ем устойчивости схемы.

Рассмотренные ниже примеры подтверждают сложность проблемы аналитического поиска общего решения (9.5) для событийно управляемых гибридных систем. Анализируемые примеры двух резер вуаров (two tanks) сформулированы достаточно недавно (1996, 2005) в разных постановках. Данные задачи являются, с одной стороны, Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ типичными представителями ГС, с другой – полновесной иллюстраци ей использования теории гибридных систем в практических задачах.

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

Событийно-непрерывная ГС Система из двух сообщающихся баков приведена на рис. 9.10.

Сигнал контроллера Input v input Winput Сигнал контроллера h v Сигнал контроллеру W H h2 Сигнал контроллера Output W2 v out Рис. 9.10. Система двух баков Трубопроводы управляются клапанами Winput, W1 и W2. Подача жидкости в систему контролируется через клапан Winput, который от крывается мгновенно со скоростью входного потока:

0, t 0, v input = 400, t 0.

Клапаны W1 и W2 являются инерционными устройствами. От момента начала открытия (закрытия) до полного открытия (закрытия) требуется 80 с. Их состояние контролируется задвижкой, меняющей свое поло 9.8. Инструментально-ориентированный анализ жение от значения P = 80 (полное закрытие) до P = 0 (полное откры тие). Управление клапанами Winput, W1 и W2 осуществляется контрол лером, входящим в описываемую систему. Если через A1 и A2 обозна чить площади оснований баков, то система уравнений для уровней жидкости в баках h1 и h2 будет иметь вид ( ) h1 = A1 1 v input v12, (9.10) h2 = A2 1 ( v12 v out ).

Скорость протекания жидкости между баками v12 зависит от уровней, и в соответствии с законом Торричелли K1 ( p1 ) h1 h2 + H, h2 H, v12 = (9.11) K1 ( p1 ) h1, h2 H.

Скорость вытекания жидкости из системы v out зависит от уровня во втором баке h2 и положения задвижки p2 на клапане W2. Если кон троллер обнаруживает, что уровень жидкости во втором баке опустил ся ниже значения Lminus, поступает команда закрыть выходной клапан Wout. Если жидкость во втором баке превышает уровень Lplus, выда ется команда открыть выходной клапан. Таким образом, имеем K 2 ( p2 ) h2, h2 Lplus, v out = (9.12) h2 Lminus.

0, Индивидуальные свойства клапанов определяются функциями 6 1,85 104 e 610 p1, 0 p1 80, K1 ( p1 ) = 0, p1 = 80, (9.13) 6 2, 26 104 e5,710 p2, 0 p2 80, K 2 ( p2 ) = 0, p2 = 80.

Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ В исходном состоянии все клапаны закрыты и оба бака пусты.

В начальный момент ( t = 0 ) контроллер посылает сигнал входному клапану Winput, который мгновенно открывается, и в течение времени Time1 [с] наполняется только первый бак. По истечении времени Time1 контроллер посылает команду открыть W1 и жидкость начинает поступать во второй бак. Второе состояние сохраняется на протяжении Time2 [с]. По истечении времени Time2 открывается W2. Если кон троллер обнаруживает, что уровень жидкости во втором баке опустил ся ниже значения Lminus, поступает команда закрыть выходной клапан Wout, если жидкость во втором баке превышает уровень Lplus – выда ется команда открыть выходной клапан.

В постановку задачи входит анализ системы с параметрами Time1 = 9, Time2 = 20, Lplus = 0,94 и Lminus = 0,3. Утверждается, что время выдачи первой команды контроллера на закрытие вентиля Wout равно 303,13 с ( h2 = Lminus = 0,3 м), а десятой – 2669,9 с. На рис. 9. приведены результаты моделирования событийно-управляемой гиб ридной системы. Представленный пример демонстрирует применение инструментальных средств компьютерного анализа гибридных систем.

Рис. 9.11. Результаты моделирования 9.8. Инструментально-ориентированный анализ В частности, изучение сфокусировано на непрерывных процессах, управляемых логическими контроллерами, как классе ГС, широко рас пространенных в технологических процессах промышленных пред приятий. Эта тестовая задача отличается от других примеров гибрид ных систем в литературе (например, контроллер «поезд-ворота» или проблема парового бойлера) по нескольким причинам.

1. Данный пример – уже существующая установка, а не идеализи рованная модель. Это становится наглядно на моделях динамики кла панов, которые были определены экспериментально в университете Дортмунда. Поэтому нереальные упрощения задачи исключены.

2. Инструментальное средство должно правильно определить со бытийно-управляемое переключение (event-triggered switching). Хотя установка выглядит очень простой и физика процесса проста для по нимания, непрерывная динамика далека от тривиальности. В результа те не очевидно с первого взгляда, как может быть создана подходящая абстракция, если будут применяться инструментальные средства, ос нованные на моделях с простой динамикой (например, HyTech с ли нейным гибридным автоматом).

3. Постановка проблемы соответствует повседневной ситуации из сферы автоматизации процессов: программа должна быть проверена в отношении соответствия требованиям неформального естественного языка, основанного на знаниях семантики процессов.

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

Событийно-управляемая ГС Рассмотрим еще одну систему двух баков (рис. 9.12), которая была предложена для иллюстрации так называемой FDI-методики диагно стики событий в гибридной системе. Система состоит из двух баков T и T2, площади которых равны (S1 = S2 ). Насос P поставляет жидкий поток Q p в бак T1. Четыре переключаемых клапана v1, v 2, v 3 и v позволяют управлять потоками Q1, Q2, Q3 и Q4. Имеются три датчи ка – два датчика потоков Q p, Q1 и один датчик для определения уров ня жидкости в баке T2. Насос управляется таким образом, чтобы под держивался уровень жидкости h2 между двумя заданными уровнями:

Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ в н h2 = 0,1 м и h2 = 0,2 м. Закон, согласно которому осуществляются переключения, определяет два потока Q p = 0 и Q p = Q0 = 0,001 м3/с.

Qp h v н h Q h v2 в h Q v v Q Q Рис. 9.12. Система двух баков Потоки Q1, Q2, Q3 и Q4 заданы в соответствии с законом Торри челли:

Q1 = A1 2qh1, Q2 = A2 2q h1 h2, (9.14) Q3 = A3 2q h1, Q4 = A4 2qh2.

Пусть v1, v 2 и v 3 всегда открыты.

pr Двум действиям клапана v 4 (открытие и закрытие) соответствуют события e1 и c1 c pr e2. В зависимости от уровня жидкости в e2 e1 e2 e баках система может находиться в одном из локальных состояний с поведениями pr c4 c ci, i = 1, 4, как показано на рис. 9.13.

Непрерывные поведения ci, i = 1, 4, pr системы с учетом динамики потоков оп Рис. 9.13. Карта поведения ределяются следующими выражениями:

9.8. Инструментально-ориентированный анализ c1: h1 0,5 ;

h2 0,5 ;

c2: h1 0,5 ;

h2 0,5 ;

( ) ( ) dh1 1 dh1 = Q p Q1 Q2 ;

= Q p Q1 Q2 Q3 ;

dt S1 dt S dh2 dh 1 ( Q2 + Q3 ) ;

( Q2 ) ;

= = dt S2 dt S c3: h1 0,5 ;

h2 0,5 ;

c4: h1 0,5 ;

h2 0,5 ;

(9.15) ( ) ( ) dh1 1 dh1 = Q p Q1 Q2 ;

= Q p Q1 Q2 ;

dt S1 dt S ( ) dh2 1 dh2 ( Q2 Q4 ).

= Q p Q1 Q2 Q3 ;

= dt S1 dt S Условия переходов из одного локального состояния в другое управля ются событиями pr1: g1 (h1 ) = h1 0,5 0, pr2: g 2 (h2 ) = h2 0,5 0, e1: g3 (t ) = t 240, e2: g3 (t ) = t 380.

Результаты машинного анализа непрерывного и дискретного поведе ния событийно-управляемой ГС двух резервуаров приведены на рис. 9.14.

Рис. 9.14. Временная диаграмма и динамика дискретных переходов Г л а в а 9. ГИБРИДНЫЕ СИСТЕМЫ Диагностика состоит в следующем. Если, например, труба P3 полно стью засорена, динамика изменения состояния будет соответствовать уравнениям состояния 1, даже если уровень h1 превысит 0,5. Это от ( ) 2 сутствие перехода определяется благодаря количеству 0,5Q1 qA1, которое превысит 0,5 м.

Приведенные элементарные и промышленные примеры гибридных систем и их математическое представление показывают, что введен ные формальные определения ГС адекватно описывают гибридные модели обозначенного класса систем. Для рассмотренных примеров характерны события, которые приводят к разрывам в первых произ водных решения. Наиболее сложными и часто встречающими в прак тических задачах являются односторонние события, для которых кро ме ограничений на событийную функцию имеет место условие неотрицательности вектора фазовых переменных. Эти и другие осо бенности гибридных систем ограничивают использование сущест вующих методов обнаружения событий. В дальнейшем остановимся на этих вопросах более подробно.

Г л а в а КОРРЕКТНОЕ ОБНАРУЖЕНИЕ ДИСКРЕТНЫХ СОБЫТИЙ О бнаружение событий – ключевая задача анализа гибридных систем. Выбор шага интегрирования по критериям точности и устойчивости может оказаться несогласованным с дискретными мо ментами смены режимов ГС. Особенно опасны в этом смысле участки установления. Здесь вследствие медленно изменяющегося решения шаг интегрирования выбирается достаточно большим, что уменьшает вероятность своевременного обнаружения момента смены локального поведения гибридной системы.

10.1. ОБЛАСТИ НЕОПРЕДЕЛЕННОСТИ ГИБРИДНОЙ МОДЕЛИ Представим плоский двухшарнирный манипулятор (рис. 10.1) с кинематическими уравнениями для угловых скоростей 1, 2 и соот ветствующими углами поворота рукавов 1, 2 в виде 1 = 1, = 2.

Координаты положения ( x, y ) для ко нечной точки задаются контроллеру проек тировщиком более высокого уровня. Отме тим, что для проектировщика высокого уровня не требуется знание специфики ма нипулятора и координат ( x, y ) точек, кото рые находятся вне множества достигаемых Рис. 10.1. Манипулятор позиций манипулятора, когда правая часть с двумя степенями свободы исходного дифференциального уравнения становится неопределенной. В этом случае, задав длину рукавов l1 и l2, l1 l2, получим событийную функцию:

( )( ) x 2 + y 2 (l1 + l2 ) x 2 + y 2 (l1 l2 ), g ( x, y ) = Г л а в а 10. КОРРЕКТНОЕ ОБНАРУЖЕНИЕ ДИСКРЕТНЫХ СОБЫТИЙ где x = l1 cos 1 + l2 cos ( 1 + 2 ), y = l1 sin 1 + l2 sin ( 1 + 2 ).

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

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

10.2. ПРОБЛЕМА КОРРЕКТНОГО ОБНАРУЖЕНИЯ ДИСКРЕТНЫХ СОБЫТИЙ На рис. 10.2 изображена диаграмма состояний некоторой гибрид ной системы. В начальный момент времени t0 система находится в состоянии c1, и непрерывное поведение системы соответствует диф ференциальным уравнениям x = f1 ( x) с начальными условиями x0 = x(t0 ). В состоянии c1 система находится, пока предикат pr : g ( x(t ) ) 0 является ложным. Как только условие g ( x(t ) ) 0 ста новится истинным, происходит переход из состояния c1 в c2. Система мгновенно переходит в другой режим, где функционирует в соответст вии с уравнением x = f 2 ( x) и теми начальными условиями, к которым система пришла в момент переключения.

g(t, x) с с x = f2(x) x = f1(x) g(t, x) Рис. 10.2. Гибридная система с двумя локальными состояниями 10.2. Проблема корректного обнаружения дискретных событий Возникает вопрос корректного обнаружения дискретных перехо дов, т. е. в процессе моделирования x = f1 ( x) на интервале [t0, t*] при g ( x0 ) 0 требуется вычислить момент t *, соответствующий началу выполнения условия g ( x(t ) ) 0.

В таких симуляторах фаза детекции заключается в проверке истин ности выражения g ( x(t0 ) ) 0. Если условие не выполняется, происхо дит численное интегрирование дифференциальных уравнений на один шаг вперед по времени для t1 = t0 + h и с последующей проверкой не равенства g ( x(t1 ) ) 0. Эта процедура повторяется до тех пор, пока не будет выполнено неравенство g ( x(tk ) ) 0. В этой точке tk считается, что событие произошло на полуинтервале (tk 1, tk ]. Отметим, что ве личина шага h выбирается без учета динамики поведения событийных функций. Некоторые системы после этого активизируют фазу локали зации для более точного определения результата, а некоторые просто принимают, что событие произошло в момент tk. Фаза локализации обычно основана на методе деления отрезка пополам. Возможны си туации, когда классический алгоритм склонен к сбоям. Примеры таких ситуаций приведены на рис. 10.3.

а б Рис. 10.3. Ситуации, в которых стандартные методы могут не обнаружить смену состояния Первый вариант изображен на рис. 10.3, а, где траектория полно стью пересекается так, что событийная функция имеет несколько кор ней на полуинтервале t* (tk 1, tk ]. Похожая ситуация возникает в тех Г л а в а 10. КОРРЕКТНОЕ ОБНАРУЖЕНИЕ ДИСКРЕТНЫХ СОБЫТИЙ случаях (рис. 10.3, б), когда множество значений событийной функции тонкое или имеет острые углы. Два этих случая практически эквива лентны, и в обеих ситуациях большинство стандартных методов могут дать сбой.

Рассмотрим ГС, в которой в зависимости от значения переменной y выбираются разные режимы. При y 0 выбирается вариант dx / dt = a1x + b1, y = x, a1 = 1, b1 = 1, x(0) = 1, а при y 0 – следующий:

dx / dt = a2 x + b2, y = x, a1 = 1, b2 = 1.

Здесь хорошо видно, что точкой переключения является значение y = 0 и продолжение решения невозможно ни слева, ни справа. Иллю страция этого случая показана на рис. 10.4. Так как правая часть ОДУ не может быть вычислена в новой точке, метод деления пополам нель зя применять для более точной локализации корня. В этой ситуации большинство стандартных методов локализации событий также не срабатывает.

P = false с1 с x = f1(x) x = f2(x) Рис. 10.4. Ситуация, когда правая часть ОДУ определена не на всем фазовом пространстве В «мультиагентной» модели ГС поиск точки переключения мето дом деления может приводить к огромным вычислительным затратам.

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

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

Эта проблема становится еще более актуальной для режимов ГС с повышенной жесткостью. Наиболее опасна для моделирования ситуа ция, когда переходный участок решения лежит вблизи границы облас ти неопределенности и якобиан событийной функции резко возрастает вблизи границы. Это может привести к «проскакиванию» точки пере ключения с большей вероятностью, чем в гладких системах.

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

10.3. ЛИНЕАРИЗАЦИЯ И МЕТОД УСТАНОВЛЕНИЯ В ЛОКАЛИЗАЦИИ СОБЫТИЙ В настоящее время существуют два конструктивных подхода для построения алгоритма, который гарантирует правильное решение за дачи Коши с ограничениями вида (9.3). Первый алгоритм основан на идее линеаризации, по аналогии с методом линеаризации в системах автоматического управления. Второй алгоритм основан на методе ус тановления и связан с поиском корней событийной функции g ( x) = на решении классической задачи Коши без ограничений.


Линеаризация событийной функции В отличие от (9.3) рассмотрим автономный режим ГС с ограниче ниями на событийную функцию dx = f ( x), dt g ( x) 0, x(0) = x0.

Г л а в а 10. КОРРЕКТНОЕ ОБНАРУЖЕНИЕ ДИСКРЕТНЫХ СОБЫТИЙ Величину шага h(k ) = dt / dk будем рассматривать как вход, а собы тийную функцию как выходное уравнение. Предположим, что вместо принадлежности множеству положительных целых чисел отсчет шагов k принадлежит непрерывному отрезку значений k [0, ). Далее до пустим, что tk – непрерывная функция вещественной переменной k, определенная как t (k ). Теперь логично записать x ( t ( k ) ) и g ( x ( t (k ) ) ).

Аналогично дискретному варианту находим «размер шага»

h(k ) = dt / dk, который является входной величиной. Динамика собы тийной функции тогда будет определяться выражением dg g dx dt =.

dk x dt dk После подстановки dx / dt = f ( x ) получим dg g = f ( x ) h( k ), (10.1) dk x где () обозначает скалярное произведение векторов. Необходимо вы брать величину шага h( k ) таким образом, чтобы g ( x) 0 при k.

Очевидно, это условие выполняется, если в (10.1) положить g h(k ) = g ( x) f ( x), (10.2) x где 0. Действительно, подстановкой (10.2) в (10.1) получим урав нение dg / dk = g, решение которого имеет вид g (k ) = g (0) exp(k ), (10.3) где g (0) – начальное значение событийной функции. Из полученного решения следует, что при k значение событийной функции асим птотически приближается к границе поверхности g ( x) = 0, никогда не пересекая ее. Метод линеаризации событийной функции иллюстриру ется на рис. 10.5. Здесь граница события достигается с заданной точ ностью, 0.

10.3. Линеаризация и метод установления в локализации событий 1 Граница события Граница события 0 –1 – Событийная функция, g Событийная функция, g –3 – g(x, t) g(k) –5 – –7 – –9 – 5 15 0 10 197 199 Время t Количество шагов k а б Рис. 10.5. Исходная событийная функция (а) и репараметризиро ванная алгоритмом линеаризации событийная функция (б) Из (10.2) следует, что на шаг интегрирования существенно влияет параметр. В непрерывном варианте линеаризации на параметр наложено слабое ограничение 0, которое при численной реализа ции нельзя использовать для обнаружения событий с любым видом событийных функций. Поэтому ниже будут наложены более жесткие ограничения на параметр.

Метод установления в обнаружении событий Пусть требуется найти корень функции g ( x) = 0 на решении клас сической неавтономной задачи Коши вида dx = f ( x), x(0) = x0.

dt Как и в предыдущем случае, переход от независимой переменной t к новой независимой переменной k приводит к системе dx = f ( x) S, x(0) = x0, dk dg g x dt g f ( x) r, g (0) = g ( x0 ), = = dk x t dk x dt = S, t (0) = 0.

dk Выбор функции S в виде S = g ( x)[g / x] Г л а в а 10. КОРРЕКТНОЕ ОБНАРУЖЕНИЕ ДИСКРЕТНЫХ СОБЫТИЙ приводит к уравнениям g dx g ( x) f ( x), x(0) = x0, = x dk dg = g ( x) f ( x), g (0) = g ( x0 ), dk g dt g ( x), t (0) = 0, = x dk которые обеспечивают устойчивое решение системы дифференциаль ных уравнений. Однако полученное методом линеаризации решение в общем случае не совпадает с полученным методом установления, S h( k ).

Перейдем к численному анализу выбора очередного шага hn +1, ко торый обеспечивает асимптотическое приближение к границе собы тийной функции.

10.4. ОБЕСПЕЧЕНИЕ АСИМПТОТИЧЕСКОГО ПРИБЛИЖЕНИЯ К ПОГРАНИЧНОЙ ПОВЕРХНОСТИ В ЯВНЫХ РАЗНОСТНЫХ СХЕМАХ Решение задачи Коши (2.1) будем искать явными методами, кото рые в векторной форме записываются в виде yn+1 = yn + hn +1n, n = 0,1, 2,... (10.4) с заданными начальными условиями y0. Здесь hn +1 – шаг интегриро вания, n – заданная гладкая N-мерная вектор-функция, зависящая от правой части задачи (2.1). Точное решение y (tn +1 ) в окрестности точки y (tn ) имеет вид y (tn+1 ) = y (tn ) + hn +1 y (tn ) + O(h 2 ).

Пусть точное y (tn +1 ) и приближенное yn +1 решения совпадают. Пре небрегая остаточным членом, с учетом (10.4) получим y (tn ) n, n = 0,1, 2,.... (10.5) 10.4. Обеспечение асимптотического приближения к пограничной поверхности Рассмотрим одностороннюю гибридную систему, записанную в виде задачи Коши с ограничениями, т. е.

y = f ( y ), y (t0 ) = y0, g ( y, t ) 0, (10.6) где g ( y, t ) – событийная функция или нелинейный предохранитель.

Поскольку многие модели g ( y, t ), представляющие интерес, линейны, будем рассматривать их как наиболее важный класс событийных функций. Любой нелинейный предохранитель можно привести к ли нейному виду добавлением дополнительной фазовой переменной x = g ( y, t ). В результате задачу (10.6) с учетом (9.2) можно переписать в виде y = f ( y), g g x = f ( y) +, x 0.

y t Особое внимание следует обратить на выбор метода интегрирова ния. Полностью неявный метод использовать нельзя, потому что он требует вычисления f ( y ) в потенциально опасной области, где модель не определена. Поэтому здесь будем использовать явные мето ды. В этом случае событийная динамика описывается соотношением g n+1 = g ( yn + hn +1n, tn + hn +1 ).

Разлагая g n +1 в ряд Тейлора и учитывая линейность g ( y, t ) и (10.5), имеем g g g n+1 = g n + hn +1 n n + n. (10.7) y t В итоге получили зависимость g n +1 от прогнозируемого шага hn +1.

Теорема 10.1. Выбор шага по формуле g g hn +1 = ( 1) g n n n + n, (0, 1), (10.8) y t обеспечивает поведение событийной динамики как устойчивой линей ной системы, которая приближается к поверхности g ( y, t ) = 0. Кро ме того, если g ( y0, t0 ) 0, то g ( yn, tn ) 0 для всех n.

Г л а в а 10. КОРРЕКТНОЕ ОБНАРУЖЕНИЕ ДИСКРЕТНЫХ СОБЫТИЙ Доказательство. Подставляя (10.8) в (10.7), имеем g n+1 = g n, n = 0, 1, 2,.... (10.9) Преобразовав рекуррентно (10.9), получим g n +1 = n+1g0. Поскольку 1, имеет место g n 0 при n. Кроме того, из условия следует, что функция g n не меняет знак. Следовательно, при g0 будет выполняться g n 0 для всех n. Тогда событийная функция ни когда не пересечет потенциально опасную область g ( yn, tn ) = 0, что завершает доказательство теоремы.

Алгоритм обнаружения с одношаговым методом второго порядка точности Рассмотрим алгоритм выбора шага интегрирования с учетом дина мики событийной функции на примере рассмотренного в разделе 5. двухстадийного метода второго порядка точности вида yn+1 = yn + p1k1 + p2 k2, k1 = hf ( yn ), k2 = hf ( yn + k1 ).

Пусть определена величина следующего шага по точности и устойчи вости hn +1 = max hn, min(h ac, h st ), p где hn – последний успешный шаг интегрирования, h ac – шаг по точ ности, h st – шаг по устойчивости. Новый шаг h ac по точности опре деляется по формуле h ac = q1hn, где q1 является решением уравнения q1 || k2 k1 || =.

Шаг h st по устойчивости задается выражением h st = q2 hn, в котором q2 определяется из соотношения q2 vn,2 = 2, где vn,2 – оценка макси мального собственного числа матрицы Якоби. Тогда управление ша гом интегрирования с учетом точности, устойчивости и динамики со бытийной функции можно выполнить по следующему алгоритму.

Шаг 1. Вычисляется f n = f ( yn, tn ).

10.4. Обеспечение асимптотического приближения к пограничной поверхности Шаг 2. Вычисляются g n = g ( yn, tn ), g n y = g ( yn, tn ) y, g n t = g ( yn, tn ) t.

g Шаг 3. Вычисляется шаг hn +1 по формуле (10.8), причем n = f n.

) ( g p Шаг 4. Вычисляется новый шаг hn +1 = min hn +1, hn +1.

Шаг 5. Выполняется следующий шаг интегрирования.

В приведенном алгоритме выбор шага интегрирования осуществ ляется с учетом точности, устойчивости и динамики событийной функции.

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

y = v y, v = g, y y = H, vy = 0, if ( y 0) and (v y 0) then v y = v y.

Длительностью фазы отскока при достаточно больших значениях ско рости можно пренебречь по сравнению с длительностью фазы полета.

Можно также рассматривать отскок как мгновенное дискретное дейст вие, результатом которого является изменение знака вертикальной со ставляющей скорости на противоположный. Момент отскока можно определить, если следить за выполнением неравенства y 0.

В начальный момент t0 = 0 имеет место y0 = 10, v x 0 = 0, v y 0 = 0, а мячик находится в состоянии «полета» и движется согласно приве денной системе уравнений. В некоторый момент времени t *, который определяется логическим предикатом y 0, мячик мгновенно меняет направление движения (v y = v y ), и полет (взлет) продолжается, но уже с новыми начальными условиями. Результаты расчета без алго ритма обнаружения событий приведены на рис. 10.6. При масштабиро вании, можно заметить, что смена знака скорости происходит не точно в тот момент, когда мячик достигает поверхности отскока (рис. 10.7).

Г л а в а 10. КОРРЕКТНОЕ ОБНАРУЖЕНИЕ ДИСКРЕТНЫХ СОБЫТИЙ Рис. 10.6. Зависимость высоты и скорости;

v (t ) = v y (t ) Рис. 10.7. Увеличенный график в момент переключения без алгоритма обнаружения событий Рис. 10.8. Увеличенный график в момент переключения с алгоритмом замедления 10.5. Метод Адамса в обнаружении событий Второй вычислительный эксперимент проводился с алгоритмом обнаружения событий. На рис. 10.8 видны качественные изменения численного решения в момент переключения, когда смена знака скоро сти происходит точно в момент отскока мячика, а решение y (t ) экспо ненциально приближается к границе режима.


10.5. МЕТОД АДАМСА В ОБНАРУЖЕНИИ СОБЫТИЙ Рассмотрим двухшаговый метод Адамса с пока неопределенными коэффициентами 1 и 2, т. е.

yn+1 = yn + hn +1 ( 1 f n + 2 f n 1 ). (10.10) Разложим в окрестности tn приближенное решение в ряд Тейлора ( ) 12 yn+1 = yn + hn +1 f n + hn+1 f n f n + O hn+1, где f n = f ( yn ), а f n = f ( yn ) / y – матрица Якоби от векторной функции f ( y ) в точке yn. Аналогично для значения yn1 запишем () 12 yn1 = yn hn +1 f n + hn +1 f n f n + O hn.

Тогда соответствующие ряды для f ( yn+1 ) и f ( yn1 ) имеют вид () f ( yn+1 ) = f n + hn +1 f n f n + O hn +1, f ( yn1 ) = f n hn +1 f n f n + O ( hn +1 ).

Подставляя f ( yn1 ) в (10.10), получим yn+1 = yn + hn +1 ( 1 + 2 ) f n 2 hn hn +1 f n f n +…, где элементарные дифференциалы f n и f n f n вычислены на прибли женном решении yn. Точное решение y (tn +1 ) в виде ряда Тейлора имеет вид ( ) 12 y (tn+1 ) = y (tn ) + hn+1 f + hn +1 f + O hn +1, f Г л а в а 10. КОРРЕКТНОЕ ОБНАРУЖЕНИЕ ДИСКРЕТНЫХ СОБЫТИЙ где элементарные дифференциалы f и f f вычислены на точном ре шении y (tn ). При условии yn = y (tn ) локальную ошибку n схемы (10.10) можно представить в виде n = y ( tn +1 ) yn +1. Учитывая пред ставления точного и приближенного решений, ее еще можно записать следующим образом:

n = Bhn +1 f n + Cf n f n, где B = 1 1 2, C = 0,5hn +1 + 2 hn hn+1.

Тогда, полагая C = B = 0, получим 2 = 0,5hn+1 / hn и 1 = 1 + 0,5hn+1 / hn.

С учетом (10.4) и (10.10) для двухшаговой схемы Адамса имеет место n = 1 f n + 2 f n1. Подставим n, 1 и 2 в (10.8), получим полином относительно текущего шага интегрирования Pn (hn +1 ) = a2 hn +1 + a1hn+1 + a0, (10.11) где 1 g a2 = y ( f n f n 1 ), 2hn g g a1 = fn +, a0 = (1 ) g n.

y t Алгоритм управления шагом в методе Адамса отличается от приве денного выше алгоритма для схемы Рунге–Кутты. Очередной шаг hn + определяется вещественными корнями полинома (10.11), причем если корни разные (hn +1,1 hn+1,2 ), то шаг выбирается как min(hn+1,1, hn +1,2 ).

Приведем алгоритм на основе метода Адамса. Пусть решение yn в точке tn вычислено с шагом hn и задан минимальный шаг h min и зна чение herr.

Шаг 1. Вычислить f n = f ( yn, tn ).

Шаг 2. Вычислить g n = g ( yn, tn ), g n y = g ( yn, tn ) y, 10.5. Метод Адамса в обнаружении событий g n t = g ( yn, tn ) t.

Шаг 3. Найти положительные вещественные корни r полинома P (hn +1 ) из (10.11). Если r = 0, то h =, иначе h = min(r ).

Шаг 4. Выбрать величину шага hn +1 = max[min( h, herr ), h min ].

Шаг 5. Выполнить следующий шаг интегрирования.

Пример. Рассмотрим движение управляемого контроллером объекта (тележки-робота) в ограниченном пространстве (коридоре), как это по казано на рис. 10.9. Кинематические уравнения движения имеют вид x cos() u y = sin() 1, u 0 где x и y – координаты движущегося объекта, u1 и u2 являются на правленной вперед скоростью и интенсивностью поворачивания, – угол поворота объекта.

Рис. 10.9. Объект, управляемый контролле ром в ограниченном пространстве Детали управления роботом и история u1 и u2 здесь опущены, но предполагается, что это обеспечивается контроллером. Задача состоит в проверке эффективности контроллера и возможности столкновения робота с преградами. Для простоты игнорируем физические размеры робота и работаем с ним, как с точкой. Поэтому предохранитель для моделирования задан системой уравнения стен вида ( y 0,5 0) (2,9 x 0).

Г л а в а 10. КОРРЕКТНОЕ ОБНАРУЖЕНИЕ ДИСКРЕТНЫХ СОБЫТИЙ На рис. 10.10 показана ситуация, когда система управления робо том не сможет избежать столкновения. Точки интегрирования вычис лены, как оказалось, за пределами «предохраняющей» области. Поэто му при моделировании не будет обнаружено столкновение робот столкнется со стенами возле поворота ( x = 2,9, y = 0,5 ).

Y x Рис. 10.10. Результат расчета без алгоритма обнаружения событий На рис. 10.11 и 10.12 иллюстрируется работа экстраполяционного алгоритма. Видно, как алгоритм интегрирования замедляется по мере приближения к пограничной поверхности.

Y x Рис. 10.11. Результат расчета с алгоритмом обна ружения событий 10.6. L-устойчивый метод в обнаружении событий Y x Рис. 10.12. Увеличенный участок замедления Рис. 10.13. Результаты расчета модели с роботом в системе CHARON Результаты можно сравнить с результатами расчета описанной мо дели в работе [261] и полученной авторами идеи применения экстра поляционного многочлена для задачи обнаружения события в системе CHARON (рис. 10.13).

10.6. L-УСТОЙЧИВЫЙ МЕТОД В ОБНАРУЖЕНИИ СОБЫТИЙ В качестве примера рассмотрим неявный метод Эйлера, который применительно к задаче (2.1) имеет вид yn+1 = yn + hf ( yn +1 ).

В этом случае событийная динамика описывается соотношением Г л а в а 10. КОРРЕКТНОЕ ОБНАРУЖЕНИЕ ДИСКРЕТНЫХ СОБЫТИЙ ) ( p p g n+1 = g yn + hn +1 f n +1, tn + hn +1, p где f n +1 = f ( yn +1 ), hn +1 – прогнозируемый шаг интегрирования. Заме тим, что f n +1 вычисляется в потенциально опасной точке. Поэтому проверку событийной динамики будем осуществлять (1,1)-методом первого порядка точности типа (7.3), т. е.

yn+1, 1 = yn + k1, Dn k1 = hf ( yn ), Dn = E ahf n, где E – единичная матрица, f n = f ( yn ) / y. Нетрудно убедиться, что ряд Тейлора для приближенного решения yn+1, 1 имеет вид () p p yn+1, 1 = yn + hn +1 f n + O hn +1.

В результате событийную функцию g n +1 можно записать так:

(( ) ) p p p g n+1 = g yn + hn +1 f n + O hn +1, tn + hn +1.

Разлагая g n +1 в окрестности точки ( yn, tn ) в ряд Тейлора и учитывая линейность g ( y, t ), имеем p g g g n+1 = g n + hn+1 n f n + n, y t где g n g ( yn, tn ) g n g ( yn, tn ) g n = g ( yn, t n ), = =,.

y y t t p В итоге получили зависимость g n +1 от прогнозируемого шага hn +1, аналогичную (10.7). Следовательно, для выбора шага можно использо вать ранее полученную формулу g g p = ( 1) g n n f n + n. (10.12) hn + y t 10.6. L-устойчивый метод в обнаружении событий Алгоритм интегрирования L -устойчивым методом с учетом про гноза шага через событийную функцию формулируется следующим образом.

Шаг 1. Вычисляются g n g ( yn, tn ) = g n = g ( yn, t n ),, y y g n g ( yn, tn ) =.

t t p Шаг 2. Вычисляется шаг hn +1 по формуле (10.12).

Шаг 3. Вычисляется величина v( jn ) по формуле v( jn ) = D1 jn (k2 k1 ), 1 jn 2.

n ac ac Шаг 4. Вычисляется шаг hn +1 по формуле hn +1 = qhn, где параметр q определяется из уравнения q 2 v( jn ) =, где требуемая точность расчетов.

Шаг 5. Вычисляется новый шаг hn +1 по формуле ( ) p ac hn +1 = min hn +1, hn +1, Шаг 6. Выполняется следующий шаг интегрирования.

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

Г л а в а 10. КОРРЕКТНОЕ ОБНАРУЖЕНИЕ ДИСКРЕТНЫХ СОБЫТИЙ x1(t) x2(t) Рис. 10.14. Система осциллирующих на пружинах грузов Системы дифференциальных уравнений для состояний «раз дельно» и «вместе» имеют различный вид. При условии ( stickness abs ( (m1 + m2 )a1 ) ) имеем dx1 dt = v1, dv1 dt = k1 (n1 x1 ) m1, a1 = k1 (n1 x1 ) m1, dx2 dt = v2, dv2 dt = k2 (n2 x2 ) m2, a2 = k2 (n2 x2 ) m2, а при условии ( x1 = x2 ) and (v1 v2 ) можно записать a1 = ( k1n1 + k2 n2 x1 (k1 + k2 ) ) (m1 + m2 ), dv1 dt = ( k1n1 + k2 n2 x1 (k1 + k2 ) ) (m1 + m2 ), dx1 dt = v1, a2 = ( k1n1 + k2 n2 x1 ( k1 + k2 ) ) ( m1 + m2 ), dv2 dt = ( k1n1 + k2 n2 x1 ( k1 + k2 ) ) (m1 + m2 ), dx2 dt = v2, d (stickiness) / dt = stickiness, где m1 и m2 – массы грузов, k1 и k2 – жесткости пружин, n1 и n2 – нейтральные координаты грузов, x1 и x2 – координаты грузов, v1 и v2 – скорости грузов, a1 и a2 – ускорения грузов, stickiness – совокупная жесткость пружин в состоянии «вместе».

Система переходит из состояния «раздельно» в состояние «вместе»

(рис. 10.15) при выполнении условия ( x1 = x2 ) and (v1 v2 ).

10.6. L-устойчивый метод в обнаружении событий ( stickness abs ( ( m + m2 ) a1 ) ) «Раздельно» «Вместе»

( x1 = x2 ) and (v1 v2 ) Рис. 10.15. Диаграмма состояний для модели двух масс Так как в состоянии «вместе» тела ведут себя как единое целое, то, со гласно закону сохранения импульса, скорости тел при переходе изме няются по формуле v1 = v2 = (m1v1 + m2v2 ) / (m1 + m2 ).

В этом состоянии массы удерживаются с экспоненциально убывающей жесткостью stickiness. Когда разность сил натяжения пружин превос ходит жесткость, система переходит в состояние «раздельно». Преди кат перехода в состояние «раздельно» имеет вид stickness abs((m1 + m2 ) a1 ).

В решении, представленном на рис. 10.16, справа хорошо видны участки, соответствующие разным состояниям. Двигаясь раздельно, массы сталкиваются и какое-то время двигаются как единое целое. Это видно по участку, где графики сливаются. Такое движение продолжа ется, пока сила натяжения пружин вновь не разделит их.

x x x2 x Time Time а б Рис. 10.16. Динамика без алгоритма обнаружения (а) и с алгоритмом (б) Г л а в а 10. КОРРЕКТНОЕ ОБНАРУЖЕНИЕ ДИСКРЕТНЫХ СОБЫТИЙ Рис. 10.17. Масштабированный участок фазовой траектории без алгоритма замедления Рис. 10.18. Масштабированный участок фазовой траектории с экстраполяционным алгоритмом замедления 10.7. Обнаружение событий инструментальными средствами а б Рис. 10.19. Масштабированный участок фазовой траектории в момент перехода из состояния «раздельно» в состояние «вме сте» без алгоритма замедления (а) и с экстраполяционным алго ритмом замедления (б) Из рис. 10.16 хорошо видно, что расчет модели без алгоритма об наружения приводит к тому, что уже через несколько смен состояний фазовые траектории качественно различаются. На рис. 10.17–10. приведены масштабированные участки графиков для лучшей иллюст рации работы экстраполяционного алгоритма замедления шага.

10.7. ОБНАРУЖЕНИЕ СОБЫТИЙ ИНСТРУМЕНТАЛЬНЫМИ СРЕДСТВАМИ В качестве тестового примера использования инструментальных средств (см. таблицу) обратимся к классической задаче двух баков. Бу дем рассматривать точки переключения при наступлении первого и десятого событий. При этом в разных системах моделирования гиб ридной системы, рассмотренных выше, события наступают в разные моменты времени. Точное время наступления первого события состав * * ляет t1 = 303,13, десятого – t10 = 2669,9. Из сравнения результатов видно, что самые плохие результаты продемонстрированы в системах SIMULINK с использованием метода Дорманда–Принса и SHIFT с применением метода Рунге–Кутты четвертого порядка.

Г л а в а 10. КОРРЕКТНОЕ ОБНАРУЖЕНИЕ ДИСКРЕТНЫХ СОБЫТИЙ Моменты обнаружения событий Шаг Точность 1-е событие (t) 10-е событие (t) MATLAB/Taylor (Рунге–Кутты четвертого порядка) 0,1 – 303,133 2670, Переменный 1е–6 303,127 2669, SIMULINK/StateFlow (Дорманда–Принса ode45) 1 – 307 0,1 – 303,5 2671, Переменный 1e–6 303,76 2670, GPROMS (Рунге–Кутты–Фельберга) Переменный 1е–4 303,21 2670, Переменный 1е–6 303,13 2669, DYMOLA (Адамса–Башфорта–Мултона) Переменный 1е–6 303,126 2669, SHIFT (Рунге–Кутты четвертого порядка) 1 – 301 0,1 – 302,3 2666, MVS Автомат 1е–6 303,127 2669, AnyLogic Автомат 1е–6 303,129 2669, ИСМА Фельберга 5п/Пер 1е–6 303,12 2670, Адамса/Переменный 1е–6 303,119 2669, L-устойчивый MK22/Переменный 1е–6 303,119 2670, Анализ модели с включенным алгоритмом обнаружения показыва ет значительно лучшие результаты в связи с учетом точности, устой чивости и событийной функции при выборе шага интегрирования.

Следует также отметить, что непрерывные режимы рассмотренной гибридной системы являются гладкими. Поэтому явный метод Фель берга переменного порядка и шага с контролем устойчивости и L-устойчивый метод МК22, эффективно настроенные для анализа же стких режимов, не показывают в данном случае своих преимуществ, хотя по точности не уступают методам с контролем только точности.

Г л а в а АДАПТИВНЫЙ МЕТОД ИССЛЕДОВАНИЯ РЕЖИМОВ ГИБРИДНЫХ СИСТЕМ ПОВЫШЕННОЙ ЖЕСТКОСТИ Д ля решения ряда задач повышенной жесткости применение явных методов, в том числе и алгоритмов интегрирования переменного порядка и шага, не представляется возможным. Вычисли тельный процесс с явными методами может потребовать гораздо больших временных затрат по сравнению с применением неявных схем. Однако явным методом можно рассчитывать переходные участ ки даже достаточно жестких задач. Разумное комбинирование явных и L -устойчивых методов по критерию устойчивости с учетом жесткости может улучшить показатели эффективного использования численных схем. Для анализа режимов ГС с односторонними событиями и режи мами повышенной жесткости предлагается применять адаптивный ал горитм с автоматическим контролем жесткости.

11.1. ОБНАРУЖЕНИЕ ЖЕСТКОСТИ Жесткость режимов гибридных систем можно обнаружить с помо щью оценки максимального собственного значения матрицы Якоби исходной задачи. Для L -устойчивых методов оценку максимального собственного числа n, max матрицы Якоби f ( yn ) y системы y = f ( y ) на n -м шаге интегрирования вычисляем через ее норму по формуле N fi ( yn ) n, max = max y j.

1i N j = Для явных методов такая оценка определяется степенным методом предложенным выше способом через ранее вычисленные стадии. Та ким образом, зная значение n, max, на каждом шаге можно контроли Г л а в а 11. АДАПТИВНЫЙ МЕТОД ИССЛЕДОВАНИЯ ГИБРИДНЫХ СИСТЕМ ровать устойчивость численной схемы с помощью неравенства h n, max D, где h – шаг интегрирования, а постоянная D ограни чивает интервал устойчивости.

Рассмотрим в связи с этим пример жесткого режима гибридной системы, заданного системой уравнений Ван-дер-Поля:

[(1 y1 y2 ) y2 y1 ], y1 = y2, y2 = eps y1 (0) = 2, y2 (0) = 0.

Степень жесткости в данном режиме задается параметром eps ( 0,1]. Как видно из результатов расчета (рис. 11.1 и 11.2), на уча стках установления метод работает вблизи границы области устойчи вости. На рис. 11.2 с логарифмическим масштабом по оси ординат приведены характеристики n, max D, где D = 3,6 – интервал устой чивости метода Фельберга шестого порядка точности, рассчитанный в разделе 5.3.

Таким образом, обнаружение жесткости сводится к контролю не равенства h n, max D. (11.1) Рис. 11.1. Первая компонента решения уравнения Ван-дер-Поля 11.2. Неявный метод с контролем жесткости Рис. 11.2. Обнаружение жесткости В практических расчетах на каждом шаге необходимо вести кон троль выполнения (11.1) с подсчетом количества нарушений данного неравенства. При превышении заданного предела производится пере ключение с явной схемы на неявный метод RADAU5. Подсчет числа нарушений (11.1) позволяет сгладить грубость оценки максимального собственного числа в явных методах.

11.2. НЕЯВНЫЙ МЕТОД С КОНТРОЛЕМ ЖЕСТКОСТИ Известно, что неявный метод RADAU5 – один из эффективных ал горитмов для решения жестких задач. Программная реализация алго ритма RADAU5 предназначена для численного решения жестких и дифференциально-алгебраических задач вида My = f (t, y ), y (t0 ) = y0, где f : R R N R N, y R N – некоторые вектор-функции, t – неза висимая переменная, M – постоянная квадратная матрица размерно сти N. Если M – сингулярная матрица, то данная задача из диффе ренциальной переходит в дифференциально-алгебраическую задачу.

В этом случае необходимо выбрать согласованные начальные значения y (t0 ) = y0 и y (t0 ) = y0. В программе RADAU5 реализован неявный метод Рунге–Кутты (трехстадийный метод Радо_IIА) с возможностью обработки ленточных матриц Якоби. В методе RADAU5 допускаются решения задач в одной из следующих форм:

B( y ) y = f (t, y ), y = f (t, y, y), B( y ) y = F (t, y, y).

Г л а в а 11. АДАПТИВНЫЙ МЕТОД ИССЛЕДОВАНИЯ ГИБРИДНЫХ СИСТЕМ Эффективность программы может быть повышена посредством ус тановки определенных параметров. Метод широко используется в со временных пакетах инструментально-ориентированного анализа ре жимов гибридных систем повышенной жесткости.

Как уже отмечалось, неявный метод на переходных (погранслой ных) участках становится неэффективным с точки зрения времени вы полнения, так как имеет соизмеримый шаг с явными методоми. В этом случае необходим алгоритм обнаружения точки перехода на явный метод на основе контроля устойчивости. Так как L -устойчивый метод в отличие от явного не имеет ограничений на величину шага по устой чивости, то применение неравенства (11.1) для контроля жесткости недопустимо. Поскольку в методе RADAU5 для расчетов применяется матрица Якоби в явном виде, то для контроля жесткости можно ис пользовать оценку нормы матрицы Якоби N fi ( yn ) n, max = f ( yn ) y = max y j.

1i N j = Данную оценку можно применять для переключения на явный ме тод. В случае выполнения неравенства (11.1) происходит переход на явную численную формулу. Нарушение данного неравенства в явном методе вызывает обратный переход.

11.3. ЯВНЫЙ МЕТОД ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Шестистадийная формула Рунге–Кутты–Фельберга задана с учетом (5.37) и (5.38) в виде yn+1 = yn + pi ki, i = k1 = hf ( yn ), k2 = hf yn + k1, k3 = hf yn + k1 + k2, 11.3. Явный метод переменного порядка и шага 1932 k4 = hf yn + k1 k2 + k3, (11.2) 2197 439 3680 k5 = hf yn + k1 8k2 + k3 k4, 216 513 8 3544 k6 = hf yn k1 + 2k2 k3 + k4 k5.

27 2565 4104 При значениях параметров 16 p1 = p51 =, p2 = p52 = 0, p3 = p53 =, 135 (11.3) 28561 9 p4 = p54 =, p5 = p55 =, p6 = p56 = 56430 50 формула (11.2) имеет пятый порядок точности и совпадает с методом Фельберга пятого порядка. В качестве оценки ошибки можно исполь зовать оценку Cn = n,5, полученную ранее в разделе 5.3 с использова нием вложенного метода и определенную по формуле 17 ( p5i p4i )ki, Сn = n,5 = 24 i = где – некоторая норма в RN, а коэффициенты p4i обеспечивают четвертый порядок точности и определены формулами 25 p41 = p42 = 0, p43 =,, 216 2197 p44 = p45 =, p46 = 0.

, 4104 Тогда для контроля точности и при выборе шага схемы (11.2) с пара метрами (11.3) можно применять неравенство Cn, а для контроля устойчивости – следующее:

Vn D, Г л а в а 11. АДАПТИВНЫЙ МЕТОД ИССЛЕДОВАНИЯ ГИБРИДНЫХ СИСТЕМ где Vn – оценка максимального собственного числа n, max матрицы Якоби, которая определятся соотношением ( 32k3 48k2 + 16k1 )i Vn = max.

( k2 k1 )i 9 1i N Постоянное число D можно взять равным 3,6, т. е. примерно равным длине интервала устойчивости схемы (11.2) с параметрами (11.3).

При использовании численной формулы (11.2) первого порядка точности с более широкой областью устойчивости применяются коэф фициенты вида p1 = 0, 41975960186956, p2 = 0, 44944365216575, p4 = 0,12199235635231 102, p3 = 0,12964196119220, (11.4) p5 = 0,66250690732054 104, p6 = 0,11118997045939 105, при которых ее область устойчивости расширена до 72 по веществен ной оси, D = 72. Метод первого порядка предполагается использовать на участке установления, где ошибки за счет неточности (11.2) невели ки. Поэтому для контроля точности вычислений можно использовать оценки локальной ошибки,1 и,1, которые можно записать сле n n дующим образом:

An =,1 = 2 4с2 k2 k1, n An =,1 = 1 2с2 hf ( yn +1 ) k1, n где для (11.2), (11.4) имеем с2 = 840 / D 2 ;



Pages:     | 1 |   ...   | 4 | 5 || 7 | 8 |   ...   | 9 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.