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

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

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


Pages:     | 1 || 3 | 4 |   ...   | 6 |

«Российская академия наук Российская ассоциация математического программирования Институт систем энергетики им. Л.А.Мелентьева СО РАН ...»

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

1 x1 = f1 (x1, x2, u), x2 = f2 (x1, x2, u), x1 (0) = x0, x2 (0) = x0, (5) 1 a u(t) b, 0 t T, x1 (tk ) = xk, x2 (tk ) = xk. (6) 1 Предположим, что оптимальное по быстродействию управление, переводящее систему из точки x0 в точку xk, имеет одну точку переключения для всех начальных условий из рассматриваемой области X 0. Требуется в фазовой плоскости (x1, x2 ) построить оптималь ную линию переключения (ЛП), состоящую из разных точек переключения управления, соответствующих различным начальным условиям. Следовательно, для построения ЛП в пространстве переменных (x1, x2 ) необходимо решить серию задач оптимального быстро действия с различными начальными условиями (x1 (0), x2 (0)). В свою очередь численное решение каждой задачи быстродействия можно свести к решению последовательности терминальных задач с различными значениями tk. Тогда при каждом фиксированном значении tk необходимо минимизировать невязку для заданного краевого условия:

(xi (tk ) xk )2 min. (7) V (tk ) = i t= 2. Методы решения задачи быстродействия.

Рассмотрим алгоритмы поиска времени быстродействия t - наименьшего tk, при ко тором x(tk ) = xk. Сначала необходимо определить интервал [tk, tk+1 ], содержащий t. Для этого при некотором пробном значении tk решается терминальная задача с применением методов из [5], [3] и, если минимум функционала (7) равен 0, то значение tk уменьшается, в противном случае увеличивается.

Алгоритм 1. Сначала рассмотрим наиболее простой алгоритм (метод дихотомии) для уточнения времени быстродействия, который часто используется при "ручном"поиске.

Допустим интервал [tk, tk+1 ], на котором функционал (7) достигает нуля, уже установ лен. Тогда для уточнения значения tk осуществляется следующая процедура:

1. Полагая tk+1 = (tk + tk1 )/2, найдем решение терминальной задачи (5) - (7).

2. Если V (tk+1 ) = 0, то полагаем tk = tk+1.

3. Если V (tk+1 ) 0, то полагаем tk1 = tk+1.

4. Если tk tk1, то возврат в п. 1.

Эта процедура стягивает интервал неопределенности со скоростью геометрической прогрессии, знаменатель которой равен 0.5. В нашем случае критерием останова про цедуры лучше взять выполнение неравенства V (tk ), где -допустимая погрешность выполнения краевого условия x(tk ) = xk.

Алгоритм 2. Теперь перейдем к рассмотрению другого алгоритма для поиска численного решения задачи быстродействия. Идея этого алгоритма состоит в сведении задачи быстро действия к задаче оптимизации некоторого параметра в управляемой системе. Отыскание оптимального значения этого параметра ведется параллельно с процессом оптимизации управления. В качестве итерационного метода оптимизации при этом используется гради ентный метод, который является обобщением метода условного градиента, а полученное в результате численное решение будет содержать как оптимальное по быстродействию управление так и наименьшее значение оптимизируемого параметра - времени быстродей ствия.

Чтобы время быстродействия стало параметром управляемой системы, необходимо вы полнить замену переменной t на переменную [0, 1]: t = p. Тогда, очевидно, dt = pd и, следовательно, правые части управляемой системы (5) следует умножить на параметр p:

dx1 dx = pf2 (x1, x2, u), x1 (0) = x0, x2 (0) = x0. (8) = pf1 (x1, x2, u), 1 d d Задача быстродействия теперь будет состоять в минимизации параметра p: p min с одновременной минимизацией функционала - невязки (7):

I(uk, pk ) = p + K (xi (tk ) xk )2 min. (9) i i= где K - весовой коэффициент (штрафной).

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

Функция Понтрягина в данной задаче принимает следующий вид H = 1 pf1 (x1, x2, u)+ 2 pf2 (x1, x2, u), где сопряженные переменные 1, 2 являются решением системы H H (10) 1 =, 2 =, x1 x 1 (1) = 2K(x1 (tk ) xk ), 2 (1) = 2K(x2 (tk ) xk ). (11) 1 Градиент функционала I(uk, pk ) теперь будет состоять из двух компонент: производ ной по управлению и производной по параметру. Формулы, по которым вычисляются эти производные, следующие:

H H k k k k Iu (u, p ) = Ip (u, p ) = (12),.

u p На управлении uk (t) на отрезке t [0, 1] вычислим градиент терминального функцио нала (9). Для этого с управлением u = uk (t) от 0 до 1 проинтегрируем систему (8) и в каж дом узле интегрирования tj [0, 1] запомним значения фазовых координат (xk (tj ), xk (tj )).

1 В последнем узле вычисляются краевые условия для сопряженной системы (11). Затем в обратном времени от момента 1 до 0 при заданном uk (tj ) и при сохраненных ранее значениях (xk (tj ), xk (tj )) интегрируем сопряженную систему (10). В каждом узле tj [0, 1] 1 вычисляется градиент Iu (uk, pk, tj ).

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

uk+1 (tj ) = uk (tj ) Iu (uk, pk, tj ), tj [0, 1].

Однако движение вдоль антиградиента может привести к нарушению ограничений на управление. Чтобы построить направление, которое при значениях параметра [0, 1] сохраняет допустимость управления uk+1 (tj ), для всех tj [0, 1], сначала найдем решение линейной задачи Iu (uk, pk, t)(u(t) uk (t))dt + Ip (uk, pk )(p pk ) min, (13) при ограничениях 1 u 1, pmin p pmax. При данных ограничениях решением этой задачи будет управление uk :

k kk u (tj ) = 1, если Iu (u, p, tj ) 0;

uk (tj ) = 1, если Iu (uk, pk, tj ) 0;

tj [0, 1], и параметр pk :

k kk p = pmax, если Ip (u, p ) 0;

pk = pmin, если Ip (uk, pk ) 0.

Значение этого параметра определяется по окончании интегрирования сопряженной системы, когда будет вычислено значение интеграла для Ip (uk, pk ). Граничные значе ния: pmin = 0(начальный момент времени), а значение pmax должно быть больше времени быстродействия.

Штрафной коэффициент K для функционала определяется точностью выполнения краевого условия.

Последовательность приближений для поиска оптимального решения определяется формулами: pk+1 = pk k (k pk ), uk+1 (tj ) = uk (tj ) + k (k (tj ) uk (tj )), tj [0, 1], k = p u 1, 2, :. Значение параметра k находится процедурой одномерного поиска. Отметим, что значение этого параметра будет общим для обоих компонент градиента и оно определя ется из условия наискорейшего убывания функционала, т.е.: k = arg min I(uk + [k u uk ], pk + )(k pk )).

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

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

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

Алгоритм 3. Метод линеаризации ограничений. На k + 1-ой итерации рассматриваемого метода решается задача минимизации заданного функционала на решениях системы (5) при линейных ограничениях. При этом для численного решения формируется задача ма тематического программирования с линейными ограничениями специальной структуры [4].

Аппроксимирующая задача строится на сетке t0, t1,..., tN, которая может быть и нерав номерной. Управление uk (tj ) и траектория xk (tj ), j = 1, N, полученные на к-ой итера ции внешнего метода оптимизации, запоминаются в узлах этой сетки, а затем в их окрест ности выполняется линеаризация ограничений.

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

Как уже было отмечено, при построении матрицы-якобиана A можно использовать технологию параллельных вычислений для расчета нескольких градиентов одновременно [5].

Тогда аппроксимирующая задача представляется в виде следующей задачи математи ческого программирования:

L(z) min, Az = b, z.

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

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

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

1) линеаризация ограничений и решение вспомогательной линейной задачи. Наиболь ших вычислительных затрат требует расчет якобиана линейной системы, поэтому здесь целесообразно применить параллельные вычисления для сокращения расчетов.

2) проверка условий оптимальности для полученного на k-ой итерации решения и, если они не выполняются, то повторение п. 1).

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

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

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

Вычислительная технология расчета ЛП состоит в применении алгоритмов 2 или для некоторых заданных значений x1 (0), xk (1) и последовательно выбираемых значений x2 (0) [0, d], где d определяет область исследуемых значений. Затем с каждым заданным начальным значением x2 (0) решалась задача быстродействия. Далее для полученного оп тимального управления определялась его точка переключения tp и находились соответ ствующие ей значения фазовых координат x1 (tp ), x2 (tp ). При последовательном изменении значений x2 (0) и нахождении соответствующих решений задачи быстродействия получа ем таблицу значений x1 (tp ), x2 (tp ). Эти значения являются координатами точек фазового пространства, через которые проходит искомая линия переключения.

Для получения численного решения задачи быстродействия с помощью алгоритма в среднем решалось около 10 терминальных задач, а применение алгоритма 2 сразу дает решение задачи быстродействия (без решения терминальных задач). Однако для верифи кации получаемых численных результатов иногда используется и алгоритм 1, как более простой и обладающий более высокой вычислительной устойчивостью. Построение одной линии переключения (ЛП) свою очередь требовало решения порядка 15 задач быстродей ствия. Таким образом, трудоемкость построения одной ЛП примерно равна трудоемкости решения 150 терминальных задач.

Список литературы [1] М.А.Дёмкин, Б.Е.Федунов, А.Д.Шараборов Траекторная защита самолета от ра кет класса "Воздух-воздух", атакующих из передней полусферы- Изв. РАН.,Теория и системы управления, 2004, N4, c. 150-157.

[2] А.И.Tятюшкин ППП КОНУС для оптимизации непрерывных управляемых систем. В кн.: Пакеты прикладных программ: Опыт использования. М.: Наука, 1989. с. 63-83.

[3] А.И.Тятюшкин Численные методы и программные средства оптимизации управля емых систем. Новосибирск: Наука,1992, 193 с.

[4] А.И.Тятюшкин Численные методы решения задач оптимального управления с огра ничениями на фазовые координаты - Изв. РАН. Теория и системы управления, 1998, N2, c. 127-133.

[5] А.И.Тятюшкин Многометодная технология для расчета оптимального управления Изв. РАН. Теория и системы управления, 2003, N3. c. 59-67.

SYNTHESIS OF OPTIMAL CONTROL IN ONE PROBLEM OF PURSUIT Yu.V. Ageeva, A.I. Tyatyushkin Institute of System Dynamics & Control Theory of SB RAS, Irkutsk e-mail: tjat@icc.ru Abstract. Technology of numerical solving of applied optimal control problems as usual is based on universal software which has well developed interface and rich arsenal of optimization methods. Such software allows to take into account a peculiarities of the problem under consideration by making of use diverse algorithms of improvement at dierent stages of iteration process.

Key words: optimization, optimal control, parallel computing, numerical methods К ЧИСЛЕННОМУ РЕШЕНИЮ ЗАДАЧ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ С ТЕРМИНАЛЬНЫМИ ОГРАНИЧЕНИЯМИ Е.В. Аксенюшкина Институт математики и экономики ИГУ, Иркутск e-mail: aksl@online.ru Аннотация. Рассматривается задача оптимального управления с терминальными ограни чениями типа равенства. Построена процедура варьирования управления, обеспечивающая одновременный спуск относительно функционала Лагранжа и штрафного функционала по методу наименьших квадратов. Цена улучшения связана с решением линейной системы для множителей Лагранжа. Указан обоснованный выбор штрафного параметра.

Ключевые слова: оптимальное управление, слабая вариация функционала.

Рассмотрим задачу оптимального управления (задача А) 0 (u) = 0 (x(t1 )) min, связанную с обыкновенной управляемой системой x(t0 ) = x x = f (x, u, t), при наличии терминальных ограничений типа равенств i (u) = i (x(t1 )) = 0, i = 1, m.

Отметим, что явные ограничения на управление u(t) в задаче отсутствуют. Доступ ным управлением в задаче А назовем всякую кусочно-непрерывную вектор-функцию u(t), t T = [t0, t1 ]. Доступное управление является допустимым, если выполнены все терминальные условия.

Введем следующий набор условий относительно данной задачи. Функции f (x, u, t), i (x), i = 0, m непрерывны по совокупности своих аргументов на R n Rr T вместе с производными по x, кроме того существуют непрерывные производные fu (x, u, t).

Для каждого функционала i определим функцию Понтрягина H i ( i, x, u, t) = i, f (x, u, t) вместе с сопряженной системой i = Hx (, x, u, t), (t1 ) = i (x(t1 )).

Пусть имеется доступное управление u(t), t T вместе с соответствующими решени ями x(t), i (t), i = 0, m фазовой и сопряженной систем.

Будем использовать обозначение Hu [t, u] = Hu ( i (t), x(t), u(t), t).

i i Построим метод улучшения управления u(t) в рамках задачи А на основе процедуры слабого варьирования.

Образуем семейство варьированных управлений по правилу t T, u (t) = u(t) + u(t), с параметром 0. Как известно [1-3], соответствующие приращения функционалов i допускают представление i (u ) i (u) = i (u) + oi (), i = 0, m, где i - слабая вариация функционала i, которая выражается по формуле i i (u) = Hu [t, u], u(t) dt.

T Пусть = (1,..., m ) - вектор множителей Лагранжа, µ 0 - штрафной параметр.

Введем нормальный функционал Лагранжа m L(u, ) = 0 (u) + i i (u) i= и штрафной функционал терминальных ограничений по принципу наименьших квадратов m 2 (u).

P (u, µ) = µ i 2 i= Образуем модифицированный функционал Лагранжа M (u,, µ) = L(u, ) + P (u, µ).

Построим слабые вариации полученных функционалов.

Для функционала L(u, ) функция Понтрягина имеет вид H L (, x, u, t) =, f (x, u, t), с сопряженной системой m = fx (x, u, t)T, (t1 ) = 0 (x(t1 )) i i (x(t1 )).

i= Пусть (t, ) - решение этой системы при x = x(t), u = u(t). Тогда m i i (t).

(t, ) = (t) + i= Слабая вариация функционала Лагранжа имеет вид L L(u, ) = Hu [t, u, ], u(t) dt, T где L L Hu [t, u, ] = Hu ((t, ), x(t), u(t), t).

С учетом решения сопряженной системы получаем итоговое выражение m 0 i L(u, ) = Hu [t, u] + i Hu [t, u], u(t) dt.

T i= Определим для функционала P (u, µ) функцию Понтрягина H P (, x, u, t) =, f (x, u, t), вместе с сопряженной системой m = fx (x, u, t)T, (t1 ) = µ i (x(t1 )) i (x(t1 )).

i= Введем обозначения i = i (x(t1 )), i = 1, m, Тогда решение сопряженной системы при x = x(t), u = u(t) примет вид m i i (t).

(t) = µ i= Отсюда получаем m H P (, x, u, t) = µ i H i ( i, x, u, t).

i= Таким образом, слабая вариация функционала P (u, µ) имеет вид m m i i i (u) = µ P (u, µ) = µ i Hu [t, u], u(t) dt.

T i= i= Слабая вариация модифицированного функционала Лагранжа с учетом полученных ва риаций имеет вид m m 0 i i M (u,, µ) = L(u, )+P (u, µ) = Hu [t, u]+ i Hu [t, u]+µ i Hu [t, u], u(t) dt.

T i=1 i= Введем обозначения m 0 i 1 u(t) = Hu [t, u] + i Hu [t, u], i= m i 2 u(t) = µ i Hu [t, u].

i= Выберем вариацию управления в следующем виде m 0 i u(t) = Hu [t, u] + (µi + i )Hu [t, u].

i= Тогда вариация модифицированного функционала Лагранжа примет вид m m 0 i 0 i M (u,, µ) = Hu [t, u] + (µi + i )Hu [t, u], Hu [t, u] + (µi + i )Hu [t, u] dt.

T i=1 i= Запишем вариации функционалов L(u, ), P (u, µ) в следующем виде L(u, ) = 1 u(t), 1 u(t) dt 1 u(t), 2 u(t) dt, T T P (u, µ) = 2 u(t), 2 u(t) dt 2 u(t), 1 u(t) dt.

T T Чтобы обеспечить знакоопределенность вариаций, потребуем выполнение условия 1 u(t), 2 u(t) dt = 0.

T В развернутой записи получаем m m i 0 i j i Hu [t, u], Hu [t, u] + i j Hu [t, u], Hu [t, u] dt = 0. (1) T i=1 i,j= Отметим, что если u(t), t T допустимое управление в задаче А, т.е. i = 0, i = 1, m, то условие (1) выполняется автоматически. В общем случае проведем преобразование в (1) m m m j 0 i j j Hu [t, u], Hu [t, u] dt + i j Hu [t, u], Hu [t, u] dt = 0.

T T j=1 i=1 j= Отсюда вытекает, что m m j 0 i j j ( Hu [t, u], Hu [t, u] dt + i Hu [t, u], Hu [t, u] dt) = 0.

T T j=1 i= Для выполнения этого условия выберем множители Лагранжа, i = 1, m, как решение i линейной системы m i j j i Hu [t, u], Hu [t, u] dt = Hu [t, u], Hu [t, u] dt, j = 1, m. (2) T T i= Отметим, что эта система впервые была получена в [1] на основе другого требования обращения в нуль вариаций i (u, ), i = 1, m. В данном случае система (2) обеспечивает условие знакоопределенности вариаций L(u, ), P (u, µ). Если выполняется условие (2), то получаем следующие выражения для вариаций введенных выше функционалов L(u, ) = ||1 u(t)||2 dt, T ||2 u(t)||2 dt, P (u, µ) = T Следовательно, в случае 1 u(t) = 0, 2 u(t) = 0 обеспечивается локальное (для малых 0) улучшение функционалов L(u, ), P (u, µ) для любого µ 0, L(u, ) L(u, ), P (u, µ) P (u, µ).

При этом автоматически происходит локальное уменьшение модифицированного функци онала Лагранжа.

Рассмотрим вопрос о выборе штрафного параметра µ в случае P (u, µ) 0. Выражение для слабой вариации функционала P (u, µ) имеет вид m P (u, µ) = µ2 i Hu [t, u]||2 dt.

i || T i= Определим меру уменьшения функционала P (u, µ) согласно следующему условию P (u + u, µ) P (u, µ) = P (u, µ) + o().

Отсюда получаем соотношение P (u, µ) = P (u, µ).

Следовательно m m µ2 i Hu [t, u]||2 dt = µ i i2.

|| T i=1 i= Таким образом, формула для параметра µ имеет вид m 1 i=1 i µ =. (3) m i || i=1 i Hu [t, u]|| T Такой выбор обеспечивает уменьшение штрафного функционала P (u, µ) на величину по рядка P (u, µ).

Таким образом, предлагаемый выбор параметров варьирования (1,..., m, µ) гаранти рует локальное улучшение двух характерных функционалов L(u, ), P (u, µ), связанных с задачей А. Далее, после -параметрического поиска процедура повторяется для нового управления, что приводит к итерационному методу приближенного решения поставленной задачи.

Список литературы [1] Энеев Т.М. О применении градиентного метода в задачах теории оптимального управления. -Космические исследования, 1966, т.4, вып.5, с. 651-669.

[2] Федоренко Р.П. Приближенное решение задач оптимального управления. - М.: Наука, 1978. - 408 с.

[3] Срочко В.А. Итерационные методы решения задач оптимального управления. - М.:

ФИЗМАТЛИТ, 2000. - 160 с.

ON NUMERICAL METHOD FOR SOLVING OPTIMAL CONTROL PROBLEM WITH TERMINAL CONDITIONS E.V. Aksenyshkina Institute for mathematics and economics of ISU, Irkutsk e-mail: aksl@online.ru Abstract. An optimal control problem with terminal conditions in equality form is considered. The proposed control varying procedure provides simultaneous descend with respect to Lagrange and penalty functionals by the least squares method. An improvement procedure consists of solving linear systems for Lagrange multipliers. The choice of penalty parameter is well justied.

Key words: optimal control, weak variation of functional.

УСОВЕРШЕНСТВОВАННЫЙ АЛГОРИТМ РАСЧЕТА ПРОЦЕССОВ С ФАЗОВЫМИ ПЕРЕХОДАМИ А.Ф. Албу, В.И. Зубов Вычислительный Центр РАН им. А.А. Дородницына, Москва e-mail: zubov@ccas.ru Аннотация. Предлагается новый итерационный алгоритм для решения задач с фазовыми пе реходами. Проводится сравнение предлагаемого алгоритма с модифицированными алгоритмами Якоби и Гаусса-Зейделя. Показывается эффективность нового алгоритма.

Ключевые слова: итерационный процесс, разностная схема, сходимость, фазовые переходы.

Задачам с фазовыми переходами посвящено большое количество работ, в которых предложены разнообразные и эффективные методы их решения. Хорошо себя зарекомен довал метод, предложенный в [1] и развитый в [2], [3]. В настоящей работе предлагается алгоритм, представляющий собой усовершенствование алгоритма [2], [3] и обладающий более высокой скоростью сходимости. Проиллюстрируем новый алгоритм на примере сле дующей задачи плавления и кристаллизации вещества (см. также [4]):

E 1 T (E) + F (r, t), (r, t) Q = {(r, t) : 0 r R, 0 t }, (1) = r(E) t r r r (2) E(r, 0) = Ein (r), 0 r R, E T (E) = [Tex T (E(R, t))], 0 t, (3) = 0, (E) r r r=0 r=R E1 CS, E E = S CS Tpl, S E E E+ = E +S, Tpl, (4) T (E) = [E +(L CL S CS )Tpl S ] 1 CL, E+ E, L kS, E E, kS + (E E )(kL kS )(E+ E ), E E E+, (5) (E) = k(E) = kL, E E+.

Здесь E(r, t) - теплосодержание вещества в точке (r, t);

, C, k - плотность вещества, его теплоемкость и коэффициент теплопроводности (индексы L и S указывают на принадлеж ность величины жидкой или твердой фазе);

- теплота плавления вещества;

Tpl - темпе ратура плавления вещества;

Ein (r) - начальное распределение функции теплосодержания;

- коэффициент теплообмена с окружающей средой;

F (r, t) - источник подводимого теп ла;

Tex - температура окружающей среды;

T (E) - температура вещества, (E) - функция теплопроводности.

Сформулированная задача решается с помощью метода сеток. Для аппроксимации краевой задачи (1)-(5) в области Q вводится неравномерная сетка = {ri, tj }, где r0 = t0 = 0, ri = ri1 + hi1, tj = tj1 + j, (i = 1,..., K;

j = 1,..., M ).

Работа выполнена при финансовой поддержке программы "Ведущие научные школы"(проект НШ 1737.2003.1) Использование неявной аппроксимации по времени и интегро-интерполяционного метода приводит нас к следующей системе конечно-разностных уравнений:

j j j j j j1 j E0 + a0 (E0 )T (E0 ) a0 (E0 )T (E1 ) = E0 + j F0, j j Eij + [ai (Eij ) + bi (Ei1 )]T (Eij ) bi (Ei1 )T (Ei1 ) j ai (Eij )T (Ei+1 ) = Eij1 + j Fij, j (1 i K 1), (6) j j j j j EK + [aK + bK (EK1 )]T (EK ) bK (EK1 )T (EK1 ) = j1 j = aK Tex + EK + j FK, (j = 1,..., M ).

Здесь введены следующие обозначения 4 j 8 j rK 4 j (2rk hK1 ) a0 =, aK =, bK =, h2 h2 (4rk hK1 ) 4rK hK1 (1 hK1 ) 0 K 4 j (2ri + hi ) 4 j (2ri hi1 ) ai =, bi =, 4ri hi (hi + hi1 ) + h3 h2 hi 4ri hi1 (hi + hi1 ) h3 + h2 hi i i1 i1 i Eij = E(ri, tj ), Fij = F (ri, tj ), (Eij ) = [(Eij + Ei+1 )/2], j (i = 1,..., K 1).

Система алгебраических уравнений (6) разбивается на M подсистем. Представим зависимость (4) температуры T от теплосодержания E в виде T (E) = µE +, где функции µ и определяются формулами 1 S C S, E E, E E E +, 0, µ(E) = 1 L CL, E+ E, 0, E E, E E E +, Tpl, (E) = ((L CL S CS )Tpl S ) 1 CL, E+ E.

L Введем также в рассмотрение (K + 1) - мерные векторы D(Ej ), L(Ej ), U(Ej ) и j, опреде j j j ляемые через компоненты (K + 1) - мерного вектора Ej = E0 E1... EK T соотношениями j j D0 (Ej ) = a0 (E0 ), DK (Ej ) = aK + bK (EK1 ), j Di (Ej ) = ai (Eij ) + bi (Ei1 ), (i = 1,..., K 1), j L0 (Ej ) = 0, Li (Ej ) = bi (Ei1 ), (i = 1,..., K), Ui (Ej ) = ai (Eij ), UK (Ej ) = 0, (i = 0,..., K 1), j j1 j i = Eij1 + j Fij, j K = aK Tex + EK + j FK, (i = 0,..., K 1).

Каждая j - я подсистема (j = 1,..., M ) системы (6) может быть записана в виде Eij + Di (Ej )T (Eij ) Li (Ej )T (Ei1 ) Ui (Ej )T (Ei+1 ) = i, j j j (7) (i = 0,..., K).

В [2] R.E. White предлагает два итерационных алгоритма для решения системы урав нений (7). Первый - модифицированный метод Якоби. Если (K + 1) - мерный вектор Vn = V0n V1n... VK T представляет собой приближение к вектору Ej, полученное на n-ой n итерации (V0 = Ej1 ), то модифицированный метод Якоби состоит в проведении итера ций, и на каждой итерации вектор Vn+1 находится с помощью следующих соотношений:

j Li (Vn )T (Vi1 )+Ui (Vn )T (Vi+1 )Di (Vn )(Vin )+i n n Vin+1=(1)Vin+ (8), (i = 0,..., K).

1 + Di (Vn )µ(Vin ) Vin+1 Vin Итерации проводятся, пока = max не станет меньше требуемой величины.

Vin i=0,...,K Второй алгоритм - модифицированный метод Гаусса-Зейделя. Здесь на каждой итера ции вектор Vn+1 определяется соотношениями j n+ n )T (Vi1 )+Ui (Vn )T (Vi+1 )Di (Vn )(Vin )+i n Li (V Vin+1=(1)Vin + (i = 0,..., K). (9), 1 + Di (Vn )µ(Vin ) В обоих алгоритмах параметр предназначен для ускорения сходимости, причем 0, E E, 1 (1 0 )(E S CS Tpl )S + 0, E E E+, (E) = 1, E+ E.

Здесь 0 - произвольный параметр, 1 0 2 и в соотношениях (8)-(9) вычисляется по параметрам на предыдущей итерации.

В процессе работы была обнаружена невысокая скорость сходимости итерационных процессов (8)-(9) и предложен способ его существенного ускорения [5]. Будем на каждом шаге итерации вектор Vn+1 определять как решение следующей системы уравнений Li (Vn )µ(Vi1 )Vi1 + [1 + Di (Vn )µ(Vin )]Vin+1 Ui (Vn )µ(Vi+1 )Vi+1 = n n+1 n n+ j = Li (Vn )(Vi1 ) Di (Vn )(Vin ) + Ui (Vn )(Vi+1 ) + i, n n (10) (i = 0,..., K).

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

Проверка эффективности нового итерационного алгоритма (10) проводилась на основе решения разных задач. Ниже приводятся результаты, полученные при решении задачи (1)-(5) (см. также [4]). Расчеты задачи выполнены при таких значениях параметров:

S = 8200, kS = 20, CS = 670, L = 8200, kL = 24, CL = 790, Tin (r) 293.15, Tpl = 1773.15, Tex = 293.15, = 1, = 1291666.615.

Источник энергии определялся соотношениями: F (r, t) = (r)f (t), где 3 1012, 0 t 7.5 103, 1, 0 r 1.25 103, f (t) = (r) = 7.5 103 t, 0, 1.25 103 r.

0, Расчеты проводились в области Q с R = 0.01 и = 0.45. Сетка по пространственной переменной и по времени выбиралась равномерной, но шаги сетки варьировались. Резуль таты расчетов представлены в таблицах 1-2. Таблица 1 соответствует параметрам сетки K = 100 и 2 = 0.5 103, таблица 2 - параметрам K = 500 и 3 = 6.25 105. В таблицах 0 - ускоряющий параметр;

метод Якоби помечен как "Я", метод Гаусса-Зейделя - как "Г-З", новый итерационный алгоритм - как "Н". Элементами таблицы являются времена работы процессора, затраченные для расчета варианта. Они определены по отношению ко времени, затраченному для расчета варианта с помощью нового итерационного алго ритма при = 105. Незаполненость клетки означает расходимость итераций. Из анализа данных можно сделать следующие выводы: метод Якоби сходится не для всех значений ускоряющего параметра 0 ;

метод Гаусса-Зейделя сходится, как правило, быстрее (для мелких сеток) метода Якоби, но не намного;

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

Табл. = 105 = 107 = 1011 = Я Г-З Я Г-З Я Г-З Я Г-З 1.0 2.10 1.86 4.15 2.94 7.97 6.07 11.41 8. 1.1 1.96 1.92 3.87 2.80 8.18 5.53 12.26 6. 1.2 2.19 2.14 3.78 3.44 11.17 6.72 17.19 9. 1.3 2.55 2.48 5.28 4.51 15.77 8.20 25.13 11. 1.4 3.43 3.42 7.92 5.73 26.73 10.91 41.49 15. 1.5 5.39 4.32 15.90 7.30 59.98 14.28 97.58 19. 1.6 5.45 10.06 19.23 26. 1.7 7.86 14.14 27.26 38. 1.8 12.49 22.24 43.06 59. 1.9 26.69 46.74 99.57 152. Н 1.00 1.21 1.61 1. Табл. = 105 = 107 = 1011 = Я Г-З Я Г-З Я Г-З Я Г-З 1.0 2.87 2.37 7.00 5.14 15.

52 11.02 22.10 15. 1.1 2.82 2.43 6.50 4.62 15.63 9.89 24.86 14. 1.2 3.05 2.72 6.60 4.69 29.54 9.52 50.24 13. 1.3 5.13 3.08 25.03 5.07 10.02 13. 1.4 3.52 5.76 10.96 15. 1.5 4.03 6.66 12.69 17. 1.6 4.77 8.02 14.54 22. 1.7 5.92 10.29 19.80 31. 1.8 7.99 14.65 32.70 49. 1.9 13.67 27.42 65.25 103. Н 1.00 1.47 1.93 2. Интересно сравнение рассмотренных алгоритмов по числу требуемых итераций. Ме тод Гаусса-Зейделя для достижения требуемой точности затрачивает меньшее число ите раций, чем метод Якоби. Новый итерационный алгоритм требует существенно меньшего числа итераций. С ростом требуемой итерационной точности число необходимых итера ций в новом итерационном алгоритме растет незначительно, в то время как в методах Якоби и Гаусса-Зейделя этот рост значителен. В таблице 3 представлено среднее необхо димое число итераций N для рассматриваемых методов. Данные приведены для расчетов, выполненных на сетке K = 500 и 3 = 6.25 105 при 0 = 1.

Преимущество нового итерационного алгоритма особенно заметно при увеличении ве личины шага по времени (таблица 4). В алгоритме Гаусса-Зейделя параметр 0 = 1.

Верхняя часть таблицы соответствует расчетам на сетке с K = 100 (шаги по времени варьируются), а нижняя - с K = 500. Все времена в таблице относительны, нормировка в верхней и нижней частях таблицы своя. Все представленные результаты относятся к случаю геометрии с осевой симметрией. Преимущество нового итерационного алгоритма было обнаружено и при рассмотрении других задач с фазовыми переходами, например, при решении задачи с плоской геометрией, которая описана в [3].

Табл. Метод Я Г-З Н 107 N = 15 N = 9 N = 1011 N = 29 N = 20 N = Табл. = 105 = 107 = 1011 = Метод K 100 Г-З 1.86 2.94 6.07 8. 100 Н 1.00 1.21 1.61 1. 100 Г-З 0.82 1.36 2.52 3. 100 Н 0.24 0.31 0.39 0. 100 Г-З 0.61 1.01 1.86 2. 100 Н 0.09 0.11 0.15 0. 500 Г-З 2.37 5.14 11.02 15. 500 Н 1.00 1.47 1.93 2. 500 Г-З 1.28 2.95 6.52 9. 500 Н 0.13 0.16 0.23 0. 500 Г-З 1.22 2.75 6.06 8. 500 Н 0.04 0.05 0.07 0. Список литературы [1] M.E. Rose А Method for Calculating Solutions of Parabolic Equations with а Free Boundary. Math. Comput., 1960, V.14, p. 249-256.

[2] R.Е. White An Enthalpy Formulation of the Stephan Problem. SIAM J. Numer. Anal., 1982, V.19, №6, p. 1129-1157.

[3] R.Е. White А Numerical Solution of the Enthalpy Formulation of the Stephan Problem.

SIAM J. Numer. Anal., 1982, V.19, №6, p. 1158-1172.

[4] А.Ф. Албу, В.А. Инякин, В.И. Зубов Оптимальное управление процессом плавления и кристаллизации вещества. Ж. выч. матем. и матем. физ., 2004, Т.44, №8, с.1364-1379.

[5] А.Ф. Албу, В.И. Зубов О модификации одной схемы для расчета процесса плавления.

Ж. выч. матем. и матем. физ., 2001, Т.41, №9, с.1434-1443.

AN IMPROVED ALGORITHM FOR SOLVING TWO-PHASE STEFAN PROBLEMS A.F. Albu, V.I. Zubov Computing Centre of Russian Academy of Sciences, Moscow e-mail: zubov@ccas.ru Abstract. A new iterative algorithm is proposed for solving problems with phase transitions formulated in terms of enthalpy. The proposed algorithm is compared with the modied Jacobi and Gauss-seidel algorithms. The eciency of the new algorithm demonstrated with performed calculations.

Key words: iterative algorithm, nite-dierent scheme, convergence, two-phase Stefan problem.

ОПТИМАЛЬНОЕ ПО БЫСТРОДЕЙСТВИЮ УПРАВЛЕНИЕ В РЕЖИМЕ РЕАЛЬНОГО ВРЕМЕНИ В.М. Александров Институт математики им. С.Л. Соболева СО РАН e-mail: vladalex@math.nsc.ru Аннотация. Предложен простой способ формирования в режиме реального времени кусочно постоянного финитного управления, переводящего линейную систему за фиксированное время из любого начального состояния в начало координат. По финитному управлению определяется оптимальное по быстродействию управление в начальный момент. Получены соотношения преобразующие последовательность финитных управлений в оптимальное по быстродействию управление. Вычисления производятся в процессе движения, сопровождения и управления системой.

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

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

В настоящей работе дано развитие метода последовательного синтеза оптимального по быстродействию управления линейными системами. Основу метода составляет моди фикация финитного управления [4].

1. Постановка задачи Пусть управляемая система описывается линейным дифференциальным уравнением x0 V.

x = A(t)x + B(t)u, x(t0 ) = x0, (1.1) Здесь x n-мерный вектор фазового состояния;

A(t) и B(t) непрерывные матрицы размеров n n и n m соответственно;

u m-мерный вектор управления, компоненты которого принадлежат классу кусочно-непрерывных функций и подчинены ограничениям |uj | Mj, Mj 0, j = 1, m. (1.2) Предполагается, что линейная система (1.1) полностью управляема, т.е.

t k rank (tk, )B( )B ( ) (tk, )d = n, (1.3) t и переводима в начало координат ограниченным управлением (1.2), т.е. x0 принадлежит области управляемости V.

Задача. Найти в режиме реального времени допустимое управление u 0 (t), переводящее си стему (1.1) из начального состояния x(t0 ) = x0 в начало координат x(tk ) = 0 за минимальное время T = tk t0.

Суть метода состоит в следующем. Формируется в режиме реального времени кусочно-постоянное финитное управление, которое способно перевести линейную систему (1.1) из заданного начального состояния x(t0 ) = x0 в начало координат x(t0 ) = 0 за фиксированное время T = t0 t0, где t произвольно. Финитное управление u0 (t) позволяет определить оптимальное по быстродействию управление u0 (t0 ) в начальный момент времени t = t0. Это управление u0 (t) сразу подается на систему (1.1), которая начинает движение, и далее происходит итерационный процесс вычисления моментов переключений оптимального управления и времени перевода системы. Через время t = t1 t0 = h, которое необходимо на вычисление нового финитного управления u1 (t), измеряются координаты фазовой точки x(t1 ). Формируется финитное управление u1 (t), которое способно перевести систему (1.1) из точки x(t1 ) в начало координат за новое фиксированное время T = t1 t1 и т.д. Доказано, что последовательность финитных управлений сходится к оптимальному по быстродействию управлению.

2. Формирование финитного управления 2.1. Определение начальных условий сопряженной системы и моментов переклю чений управления. На интервале [t0, tk ], где tk произвольно, выбираем произвольно n p (n 1) моментов переключений j, j = 1, m;

p = 1, qj ;

qj = n 1 для компонент вектора i= управления. Для линейной системы (1.1) оптимальное по быстродействию управление задается выражением u0 (t) = Mj sign [Bj (t)] (t), j = 1, m, где (t) = (t, t0 )(t0 ) реше j ние сопряженной системы = A (t). Моменты переключений и их число однозначно определяются функциями переключений [Bj (t)] (t, t0 )(t0 ), если известно начальное условие (t0 ) сопряженной системы. В заданные моменты переключений функция пе реключений равна нулю. Получаем систему из (n1) линейных алгебраических уравнений m p p [Bj (j )] (j, t0 )(t0 ) qj = n = 0, j = 1, m;

p = 1, qj ;

(2.1) j= с (n 1) неизвестными (t0 ), которая связывает (n 1) произвольно заданных моментов переключений с начальными условиями нормированной сопряженной системы. Заданием (n 1) моментов переключений однозначно задаются на произвольном интервале [t 0, tk ] все остальные моменты переключений (если они существуеют).

2.2. Вычисление весовых коэффициентов финитного управления. Финитное управле ние формируется по алгоритму n p p1 p u (t) = Nij xi (t0 ) sign [Bj (t)] (t), j = 1, m;

p = 1, rj ;

t [j, j ]. (2.2) j i= p С учетом обозначения Nij = sign [Bj (t)] (t), получаем для финитного управления следу ющее простое выражение n p p1 p u (t) j = 1, m;

p = 1, rj ;

t [j, j ].

= Nij xi (t0 ), (2.3) j i= p Весовые коэффициенты Nij находятся из решения систем линейных алгебраических урав нений p j rj m p 1 (, t0 )Bj ( )d + Ii = 0, Nij i = 1, n. (2.4) j=1 p=1 p j 2.3. Выбор опорной последовательности. Проверяем смену знаков у финитных управ ляющих воздействий up, j = 1, m и выделяем такие номера интервалов p = p, = 1, 2,..., j для которых знаки на этом и последующем интервалах совпадают. Финитное и оптималь ное управления совпадают по знакам для такой компоненты и на такой последователь ности интервалов p [pv1 + 1, pv ], для которой максимально интегральное воздействие на систему, т.е. достигается максимум выражения pv p |up |( p p |up |(j j ).

p p ) = max (2.5) j p,j, p=pv1 +1 p=p1 + Знак оптимального управления в начальный момент t0 находим путем пересчета интер валов p u0 (t0 ) = M sign [(1)(pv 1) uv ]. (2.6) Если знаки финитного управления чередуются на всех интервалах, то u0 (t0 ) = j Mj sign u (t0 ).

j Финитное управление переводит линейную систему из любого начального состояния в начало координат за фиксированное время, но могут нарушаться ограничения (1.2). Суть предлагаемого метода нахождения оптимального управления заключается в постепенном выравнивании величин финитного управления до предельных значений ±Mj, j = 1, m, с которыми формируется оптимальное по быстродействию управление. По мере выравни вания управляющих воздействий соответствующим образом изменяются и моменты пе реключений финитного управления, которые стремятся к моментам переключений опти мального управления. Однако это выравнивание должно происходить постепенно, иначе возможна расходимость вычислительного процесса из-за приближенности полученных со отношений.

p Моменты переключений j, j = 1, m, p = 1, rj 1, время перевода T = tk t0 и весовые коэффициенты Nip финитного управления не зависят от начального условия x(t0 ) = x0, j как это непосредственно видно из (2.4). Именно это свойство позволяет находить весовые коэффициенты и все моменты переключений предварительно до начала процесса управ ления. Простота реализации финитного управления (2.3) позволяет в режиме реального времени определить знаки компонент искомого оптимального управления в начальный момент t0, т.е. вычислить u0 (t0 ).

3. Итерационный процесс вычисления оптимального управления 3.1. Вариация управляющих параметров. Отклонение величины финитного управ ления от соответствующего предельного значения Mj Sj (p) = Mj sign [Bj (t)] (t), p1 p t [j, j ] при оптимальном управлении для j компоненты на p интервале составляет n p up p1 p = Mj Sj (p) j = 1, m;

p = 1, rj ;

t [j, j ].

Nij xi (t0 ), (3.1) j i= Отклонения up порождают следующее отклонение фазовых координат в конечный мо j мент tk :

p j rj m n p (tk, )Bj ( ) Mj Sj (p) (tk ) = x Nij xi (t0 ) d. (3.2) j=1 p=1 i= p j 3.2. Вариация моментов переключений управления. Вариация моментов переключений p на j, а конечного момента на tk, для кусочно-постоянного управления u(t), компонен ты которого переключаются в моменты времени j и принимают значения uj (t) = up, p j p1 p t [j, j ], вызывает следующее отклонение фазовых координат:

m rj 1 m (tk ) r (tk, j )Bj (j )[up uj ]j + p p p1 p Bj (tk )uj j tk.

x (3.3) = j j=1 p=1 j= 3.3. Уравнение баланса отклонений. Отклонение (tk ), вызванное изменением вели x чин управляющих воздействий, должно быть скомпенсировано отклонением (t k ), вы x званным изменением моментов переключений:

(tk ) + (tk ) = 0.

x x (3.4) Подставим в (3.4) выражения (3.2) и (3.3). Получим систему из n линейных алгебраиче ских уравнений. Для нахождения оптимального управления (а не просто допустимого) необходимо перейти к определению числа и расположения моментов переключений с по мощью сопряженной системы.

3.4. Связь между приращениями начальных условий нормированной сопряженной системы и приращениями моментов переключений управления. Для j компоненты вектора управления функция переключения равна нулю в каждый из (rj 1) моментов переключений. Проварьировав (t0 ) на (t0 ) в полученном уравнении (2.1) и разложив его в ряд Тейлора, найдем искомую связь:

j p p p p p p p [Bj (j )] A (j )[Bj (j )] [1 (j, t0 )] (t0 ) [Bj (j )] [1 (j, t0 )] (t0 ), = (3.5) j = 1, m;

p = 1, rj 1.

Подставив (3.5) в (3.4), получим систему из n линейных алгебраических уравнений с n неизвестными, решив которую находим (t0 ) и tk.

3.5. Вычисление оптимального управления. По формуле (3.5) вычисляем приращения p j моментов переключений управления. Проверяем выполнение условия p max[|j |, |tk |] (tk t0 ), 0 1. (3.6) j,p Если (3.6) выполняется, то находим к моменту t1 (отстоящему от t0 на время вычисления p h) уточненные значения моментов переключений j (t1 ) и конечного момента tk (t1 ):

p p p j (t1 ) = j (t0 ) + j ;

tk (t1 ) = tk (t0 ) + tk.

Вычисляем значение сопряженной системы в момент t1 :

(t1 ) = (t1, t0 )((t0 ) + (t0 )).

Находим значение фазовых координат t (t1, )B( )u0 ( )d, x(t1 ) = (t1, t0 )x(t0 ) + t которые являются начальными условиями в момент t1. Начинается новый шаг вычисления и уточнения моментов переключений управления и конечного момента.

p p p Если (3.6) не выполняется, то полагаем j (t1 ) = j (t0 ) + j ;

tk (t1 ) = tk (t0 ) + tk, где коэффициент вычисляется по формуле (tk (t0 ) t0 ) =.

p max [ |j |, |tk | ] j,p В результате максимальное отклонение на этой итерации принимается равным предельно допустимому значению, а остальные отклонения уменьшаются пропорционально в раз, и компенсируется не все отклонение (tk ), а лишь его часть. Ограничение больших x отклонений позволяет избежать расходимости итерационного процесса вычислений из-за приближенности соотношений (3.3) и (3.5). Вычисляем значение сопряженной системы в момент t (t1 ) = (t1, t0 )((t0 ) + (t0 )) и находим значение фазовых координат x(t1 ). Если в процессе движения системы (1.1) происходит переключение управления, то фазовая траектория вычисляется с учетом пе реключения управления и продолжается уточнение моментов переключений и времени процесса. Вычисление оптимального управления сводится к последовательному решению систем линейных алгебраических уравнений и задач Коши. Так постепенно осуществляет ся выравнивание величин финитного управления до предельных значений ±Mj в процессе движения, сопровождения и управления системой (1.1). Финитное управление с предель ными значениями управляющих воздействий ±Mj переводит систему в начало координат, находится с помощью сопряженной системы, удовлетворяет необходимому условию опти мальности принципу максимума Понтрягина и является оптимальным по быстродей ствию управлением.

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

Список литературы [1] Р. Габасов, Ф.М. Кириллова. Оптимальное управление в режиме реального времени.

- Вторая международная конференция по проблемам управления (17-19 июня 2003 г.).

Пленарные доклады. М.: Институт проблем управления. 2003, с. 20-47.

[2] В.М. Александров Последовательный синтез оптимального по быстродействию управления. - Журнал вычислительной математики и математической физики. 1999, т. 39, № 9, с. 1464-1478.

[3] В.М. Александров Сходимость метода последовательного синтеза оптимального по быстродействию управления. - Журнал вычислительной математики и математиче ской физики. 1999, т. 39, № 10, с. 1650-1661.

[4] В.М. Александров Приближенное решение задачи линейного быстродействия. Авто матика и телемеханика. 1998, № 12, с. 3-13.

TIME OPTIMAL CONTROL IN REAL TIME MODE V.M. Aleksandrov Sobolev Institute of Mathematics, Siberian Branch of the Russian Academy Sciences e-mail: vladalex@math.nsc.ru Abstract. It is proposed a simple method of forming in real time nite piecewise constant control which moves linear system in xed time from an initial state to the origin. The nite control is used for time optimal control to be found at initial moment. Obtained are relations transforming the succession of nite controls into time optimal control. Computations are carried out while the system monitored.

Key words: optimal control, nite control, linear system, phase trajectory, high-speed, switching moments, adjoint system, variation, iteration.

МЕТОД ЧИСЛЕННОГО РЕШЕНИЯ НЕЛИНЕЙНОЙ ЗАДАЧИ ОПТИ МАЛЬНОГО УПРАВЛЕНИЯ СПЕЦИАЛЬНОГО ВИДА В.Г. Антоник, А.Л. Ветрова Иркутский государственный университет, Иркутск e-mail: vga@math.isu.ru Аннотация. Рассматривается нелинейная задача оптимального управления, где функции, входящие в целевой функционал представлены в виде сумы двух функций, одна из которых является вогнутой. На основе квадратичной аппроксимации задачи строится метод улучшения допустимых управлений.

Ключевые слова: оптимальное управление, вогнутый функционал, метод улучшения 1. Постановка задачи Рассмотрим нелинейную задачу оптимального управления в следующей постановке F (x(t), t)dt min, (1) (u) = (x(t1 )) + T x = f (x, u, t), x(t0 ) = x0, (2) u(t) U, t T = [t0, t1 ]. (3) Предположим, что:

1) вектор-функция f (x, u, t) непрерывно-дифференцируема по x R n, непрерывна по u U и кусочно-непрерывна по t T ;

2) функции (x) и F (x, t) представимы в виде (4) (x) = 1 (x) + 2 (x), F (x, t) = F1 (x, t) + F2 (x, t), где 1 (x) и F1 (x, t) – дважды непрерывно-дифференцируемые по x Rn функции, 2 (x) и F2 (x, t) непрерывно-дифференцируемы и вогнуты по x. Кроме того, F1 (x, t) и F2 (x, t) – кусочно-непрерывны по t T ;

3) множество U Rr компактно.

Определим класс допустимых управлений V как множество кусочно-непрерывных век тор функций u(t) с условием (3).

Приведем частный случай задачи (1)-(4) – квадратично-вогнутую задачу опти мального управления:

1 (x) = x, M x, F1 (x, t) = x, N (t)x, f (x, u, t) = A(u, t)x + b(u, t).

Работа выполнена при финансовой поддержки РФФИ (проект 02-01-00243) и программы "Универси теты России"(проект ур.03.01.002).

Введем в рассмотрение функцию Понтрягина для задачи (1)-(4) H(, x, u, t) =, f (x, u, t) относительно сопряженной вектор-функции = (t), удовлетворяющей сопряженной системе = Hx (, x, u, t), (t1 ) = x (x(t1 )).

2. Аппроксимация задачи. Метод улучшения Проведем аппроксимацию задачи (1)-(4). Пусть u, v V – пара допустимых управле ний, x(t, u), x(t, w) = x(t, u)+x(t), t T – соответствующие фазовые траектории. Введем в рассмотрение векторно-матричную задачу = Hx (, x(t, u), u(t), t), P = fx (x(t, u), u(t), t)T P P fx (x(t, u), u(t), t) (5) Hxx (, x(t, u), u(t), t) + F2xx (x(t, u), t), (t1 ) = x (x(t1, u)), P (t1 ) = 1xx (x(t1, u)).

Пусть (t, u), P (t, u), t T – решение векторно-матричной задачи Коши (5). Тогда с учетом свойства вогнутости функций 2 (x), F2 (x, t) оценка приращения функционала принимает вид w (u) (6) w(t) H((t, u) + P (t, u)x(t), x(t, w), u(t), t)dt +, T где остаток определяется соотношением 2 ) = o1 ( x(t1 ) oH ( x(t) )dt T x(t), P (t, u)of ( x(t) ) dt.

T Отметим, что в квадратично-вогнутой задаче = 0, т.е. оценка (6) является нелокаль ной (точной) [2].

Далее, определим процедуру игольчатого варьирования управления u(t) с помощью произвольного управления v V и функции варьирования (t), удовлетворяющей усло виям (t) {0, 1}, (t)dt = (t1 t0 ).

T Здесь [0, 1] – параметр варьирования. Тогда оценка приращения функционала пред ставляется следующим образом w (u) 2 (u, v, ) + o(2 ), 2 (u, v, ) = (t)v(t) H(p(t, u, x(t, uv, )), x(t, uv, ), u(t), t)dt, T p(t, u, x) = (t, u) + (t, u)(x x(t, u)).


Используя этот результат, построим метод улучшения управления u(t) в задаче (1)-(4).

По заданному допустимому управлению u(t) найдем соответствующую фазовую траекто рию x(t, u), t T. Решим векторно-матричную задачу Коши (5) с целью поиска функций (t, u), P (t, u).

Далее, построим вспомогательную вектор-функцию p(t, u, x) = (t, u) + (t, u)(x x(t, u)), x Rn, t T, сформируем максимизирующее управление v (x, t) = arg max H(p(t, u, x), x, w, t), x Rn, t T.

wU Образуем функцию переключения g(x, t) = v(t,x) H(p(t, u, x), x, u(t), t).

Построим параметрическое семейство управлений u(t), g(x, t) u(x, t, ) = v (x, t), g(x, t) относительно параметра 0.

Найдем решение x (t) фазовой системы x = f (x, u(x, t, ), t), x(t0 ) = x и сформируем управление u (t) = u(x (t), t, ).

Параметр 0 найдем из условия улучшения (u ) (u).

Сделаем ряд замечаний.

1) В нелинейной квадратично-вогнутой задаче описанный метод обеспечивает нелокальное улучшение для любого значения 0.

2) Для поиска параметра можно использовать процедуру половинного деления, выбирая в качестве стартового значения величину = inf g(x(t, u), t).

tT 3) Если в задаче (1)-(3) условие (4) не выполнено, то можно воспользоваться следующей схемой:

(x) = (x) (x) + (x), F (x, t) = F (x, t) F (x, t) + F (x, t), где (x), F (x, t) – вогнутые по x Rn функции. В результате, 1 (x) = (x) (x), 2 (x) = (x), F1 (x) = F (x, t) F (x, t), F2 (x) = F (x, t).

В качестве функций (x), F (x, t) можно взять следующие (x) = x, Dx, F (x, t) = x, Q(t)x, где D 0, Q(t) 0, t T.

Список литературы [1] Срочко В. А. Итерационные методы решения задач оптимального управления. – М.:

Физматлит, 2000.

[2] Антоник В.Г., Ветрова А.Л. Итерационный метод решения вогнутой задачи опти мального управления // Вестник Бурятского университета. Серия 13: Математика и информатика. – Вып. 1. – Улан-Удэ: Изд-во Бурятского госуниверситета, 2004. – С.65-69.

NUMERICAL METHOD FOR SOME NONLINEAR OPTIMAL CONTROL PROBLEM V.G. Antonik, A.L. Vetrova Irkutsk State University, Irkutsk e-mail: vga@math.isu.ru Nonlinear optimal control problem is considered. In this problem the cost functional is presented as a sum of concave functional and a nonlinear functional. The method of control improvement is constructed.

Key words: Optimal Control, Concave Functional, Improvement Method К ВОПРОСУ ОБ ОПТИМАЛЬНОМ УПРАВЛЕНИИ НАЧАЛЬНЫМИ УСЛОВИЯМИ ЛИНЕЙНОЙ ПО СОСТОЯНИЮ КАНОНИЧЕСКОЙ ГИПЕРБОЛИЧЕСКОЙ СИСТЕМЫ А.В. Аргучинцев, Е.А.Лутковская Иркутск, Иркутский государственный университет e-mail: arguch@math.isu.ru;

elut@math.isu.runnet.ru Аннотация. В работе рассматривается задача оптимального управления линейной по состоянию канонической гиперболической системой первого порядка с начальными условиями, заданными в виде управляемой системы обыкновенных дифференциальных уравнений. Для решения задачи предлагается использовть неклассическую точную формулу приращения целевого функционала для сведения исходной задачи к значительно более простой задаче оптимального управления системой обыкновенных дифференциальных уравнений.

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

Введение Известно, что задачи оптимального управления каноническими гиперболически ми уравнениями первого порядка возникают при исследовании целого ряда химико технологических процессов [1], [2]. В работах [2], [3] получено необходимое условие опти мальности вида принципа максимума Л.С.Понтрягина для распределенных управлений, входящих в правые части системы. В [4] доказан принцип максимума для сосредоточен ных на границе области управлений, которые входят в начальные условия, задаваемые в виде систем обыкновенных дифференциальных уравнений.

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

Вместе с тем, для такого класса задач оказалось весьма эффективным использование двух нестандартных точных формул приращения целевого функционала, разработанных в [5],[6]. Такие формулы приращения справедливы для произвольных допустимых процес сов и не содержат остаточных членов. Оба варианта позволяют свести исходную задачу оптимального управления к значительно более простой задаче оптимального управления системой обыкновенных дифференциальных уравнений. При этом гиперболическую систему необходимо проинтегрировать всего лишь два раза: в начале процесса и в самом его конце. В данной статье мы остановимся на использовании только одной из формул Работа выполнена при финансовой поддержке РФФИ (проект 05-01-00187) и программы "Универси теты России"(проект УР.03.01.064) приращения [5],[6]. Второй вариант формулы приращения является симметричным и приводит к эквивалентной задаче оптимального управления системой обыкновенных дифференциальных уравнений.

1. Постановка задачи В прямоугольнике = S T, S = [s0, s1 ], T = [t0, t1 ] рассмотрим каноническую гиперболическую систему первого порядка y z = f (1) (y, z, s, t), = f (2) (y, z, s, t), (1.1) t s f (1) (y, z, s, t) = A1 (s, t)y + B1 (s, t)z + q (1) (s, t), f (2) (y, z, s, t) = A2 (s, t)y + B2 (s, t)z + q (2) (s, t), (s, t).

Здесь функции y(s, t) = (y1 (s, t), y2 (s, t),..., yn1 (s, t)), z(s, t) = (z1 (s, t), z2 (s, t),..., zn2 (s, t)) определяют состояние процесса;

A1, A2, B1, B2 – матричные функции размерности соот ветственно n1 n1, n2 n1, n1 n2, n2 n2 ;

q (1) и q (2) – соответственно n1 – мерная и n – мерная вектор–функции.

Начальные условия для системы (1.1) определяются из управляемых систем обыкно венных дифференциальных уравнений y(s, t0 ) = N (u(s), s)y(s, t0 ) + µ(u(s), s), s z(s0, t) = M (w(t), t)z(s0, t) + (w(t), t), (1.2) t y(s0, t0 ) = y 0, z(s0, t0 ) = z 0.

Здесь N и M – матричные функции размерности n1 n1 и n2 n2 соответственно;

µ(u, s) и (w, t) – соответственно n1 – мерная и n2 – мерная вектор–функции;

y 0 и z 0 – заданные вектора.

Под допустимыми управлениями будем понимать ограниченные и измеримые соответ ственно на S и T функции u(s) и w(t), удовлетворяющие почти всюду на этих отрезках ограничениям типа включения u(s) U E r1, s S;

w(t) W E r2, t T. (1.3) Целью задачи оптимального управления является минимизация функционала c(1) (s), y(s, t1 ) ds+ J(u, w) =, y(s1, t0 ) +, z(s0, t1 ) + S c(2) (t), z(s1, t) dt + b(1) (s, t), y(s, t) + b(2) (s, t), z(s, t) + ds dt, (1.4) T где c (s) и b (s, t) – n1 – мерные, c(2) (t) и b(2) (s, t) – n2 – мерные вектор-функции, и (1) (1) – заданные вектора.

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

Метод последовательных приближений, примененный в [7] для обоснования существо вания и единственности решения в широком смысле задачи Коши для системы полулиней ных гиперболических уравнений первого порядка, позволяет обосновать существование и единственность непрерывных в прямоугольнике функций y(s, t) и z(s, t), удовлетворя ющих в каждой точке (s, t) интегральным уравнениям t s f (1) (y, z, s, )d, z(s, t) = z(s0, t) + f (2) (y, z,, t)d, y(s, t) = y(s, t0 ) + t0 s где y(s, t0 ) и z(s0, t) – решения задач Коши (1.2).

В задаче оптимального управления (1.1)–(1.4) матрицы коэффициентов в (1.2) зависят от управления. Это обстоятельство является причиной того, что принцип максимума [4] здесь не будет являться достаточным условием оптимальности даже в случае линейного целевого функционала. Методы последовательных приближений [3] в данном случае не учитывают специфику задачи и имеют ту же структуру, что и в общих нелинейных задачах. К значительно более интересным результатам приводит в данной задаче при менение нестандартной точной (без остаточного члена) формулы приращения целевого функционала [5],[6].

2. Формула приращения функционала С целью уменьшения громоздкости обозначений поясним идею подхода для случая, когда управление сосредоточено только на одной из границ прямоугольника. Пусть этим управлением будет u = u(s), то есть в (1.2) функции M и не зависят от управляющих воздействий.

Рассмотрим два допустимых процесса {u;

y, z} и { = u + u;

u y = y + y, z = z + z}. Тогда система в приращениях имеет вид y = A1 (s, t)y + B1 (s, t)z, t z = A2 (s, t)y + B2 (s, t)z, (s, t) ;

(2.1) s y(s, t0 ) = N (, s)(s, t0 ) N (u, s)y(s, t0 ) + µ(u, s), uy s y(s0, t0 ) = 0, s S;

z(s0, t) = 0, t T. (2.2) Здесь µ(u, s) = µ((s), s) µ(u(s), s).

u Преобразуем правую часть формулы (2.2), используя представление N (, s)(s, t0 ) N (u, s)y(s, t0 ) = u N (u, s)(s, t0 ) + N (u, s)y(s, t0 ), uy y (2.3) запишем на основании этого представления формулу приращения целевого функционала и далее проведем традиционные при доказательстве принципа максимума преобразова ния этой формулы, отказавшись только от перехода к невозмущенным значениям y(s, t 0 ) для первого слагаемого правой части (2.3). Тогда получим точную формулу приращения целевого функционала J(u) = J() J(u) = u p(s), u N (u, s)(s, t0 ) + µ(u, s) ds, y (2.4) S в которой функция p(s) определяется из следующей сопряженной задачи:


(1) = AT (s, t) (1) AT (s, t) (2) + b(1) (s, t), 1 t (2) = B1 (s, t) (1) B2 (s, t) (2) + b(2) (s, t), (s, t) ;

T T (2.5) s (1) (s, t1 ) = c(1) (s), s S;

(2) (s1, t) = c(2) (t), t T ;

p = N T (u, s)p (1) (s, t0 ), s S, (2.6) p(s1 ) =.

Таким образом, в формулу приращения (2.4) входит решение сопряженной задачи (2.6), подсчитываемое для "старого" невозмущенного управления u = u(s), и решение задачи Коши (1.2), вычисляемое для "нового" возмущенного управления u = u + u.

3. Вариационный принцип максимума Отметим принципиальное отличие формулы (2.4) от формулы приращения, обычно используемой при доказательстве принципа максимума: формула (2.4) не содержит оста точных членов. Однако при этом система система (1.2) интегрируется на возмущенном управлении u.

Формула приращения (2.4) позволяет свести исходную задачу оптимального управле ния (1.1)–(1.4) к линейной по состоянию задаче оптимального управления системой обык новенных дифференциальных уравнений.

Теорема.Пусть u = u(s) – произвольное допустимое управление в задаче (1.1)–(1.4).

Допустимое управление u = u (s) оптимально в этой задаче в том и только том случае, когда I(u ) = min I(v), v I(v) = p(s, u), (N (v(s), s) N (u(s), s))x(s, v) + µ(v(s), s) µ(u(s), s) ds, S x = N (v, s)x + µ(v, s), (3.1) x(s0 ) = y 0, v(s) U, s S, где p(s, u) – решение задачи (2.6) при u = u(s);

x = x(s) есть n1 -мерная функция состо яния;

v = v(s) – функция управления, удовлетворяющая тем же ограничениям, что и управляющая функция исходной задачи оптимального управления.

Доказательство непосредственно следует из формулы приращения (2.4). Из этой же формулы вытекает, что оптимальное значение функционала в задаче (1.1)–(1.4) равно J(u ) = J(u) + I(u ).

Данное условие оптимальности носит вариационный, а не конечномерный характер.

Поэтому по аналогии с [3] оно названо вариационным принципом максимума.

Таким образом, для решения задачи оптимального управления (1.1)–(1.4) необходимо выполнить следующие операции.

1. Задать произвольное допустимое управление u(s). Вычислить соответствующее ему решение p(s, u) сопряженной задачи (2.6). Отметим, что для этого требуется также найти решение задачи (2.5), которое не зависит от выбора допустимого процесса.

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

Следовательно, при реализации данной схемы нужно всего лишь два раза проинте грировать уравнения с частными производными (поиск = (s, t) и при необходимости x = x(s, t, u )).

Итак, исходная задача оптимального управления в канонической гиперболической си стеме сведена к линейной по состоянию задаче оптимального управления системой обык новенных дифференциальных уравнений (3.1). Аналогичный результат может быть полу чен для управления w. Для решения исходной задачи (1.1)–(1.4) путем сведения к задаче (3.1) необходимо всего лишь два раза проинтегрировать уравнения с частными произ водными. Подчеркнем еще раз, что если бы задача (1.1)–(1.4) решалась итерационными процессами принципа максимума, то на каждой итерации приходилось бы неоднократно интегрировать гиперболическую систему (1.1). Отметим также, что для решения задачи (3.1), на данный момент уже хорошо изученной, можно использовать весь набор достаточ но эффективных методов, разработанных для задач оптимального управления системами обыкновенных дифференциальных уравнений. Весьма существенно, что эта задача опти мального управления являются линейной по состоянию, причем матрицы коэффициентов при фазовых переменных в дифференциальных уравнениях зависят от управляющих воз действий. Поэтому в данном случае эффективно применение методов [5], учитывающих указанную специфику задачи.

Список литературы [1] Островский Г.М., Волин Ю.М. Методы оптимизации сложных химико технологических систем. М.: Химия, 1970.

[2] Островский Г.М., Волин Ю.М. Моделирование сложных химико-технологических си стем. М.: Химия, 1975.

[3] Васильев О.В., Срочко В.А., Терлецкий В.А. Методы оптимизации и их приложе ния. Часть 2. Оптимальное управление. Новосибирск: Наука. Сиб. отд-ние, 1990.

[4] Васильев О.В. Принцип максимума Л.С.Понтрягина в теории оптимальных систем с распределенными параметрами // Прикладная математика. Новосибирск, Наука.

Сиб. отд-ние, 1978.

[5] Срочко В.А. Итерационные методы решения задач оптимального управления. М.:

Физматлит, 2000.

[6] Аргучинцев А.В. Оптимальное управление начально-краевыми условиями гипербо лических систем. Иркутск: Изд-во Иркутского гос. университета, 2003.

[7] Рождественский Б.Л., Яненко Н.Н. Системы квазилинейных уравнений и их прило жения к газовой динамике. М.: Наука, 1978.

ON OPTIMAL CONTROL OF INITIAL CONDITIONS OF A LINEAR WITH RESPECT TO STATE CANONICAL HYPERBOLIC SYSTEM A.V. Arguchintsev, E.A. Lutkovskaya Irkutsk State University, Irkutsk e-mail:arguch@math.isu.ru;

elut@math.isu.runnet.ru Abstract. A problem of optimal control of a linear with respect to state, canonical hyperbolic system of rst order is considered. Here the initial conditions are given by a controlled system of ordinary dierential equations. To solve this problem in the paper it is proposed to use a non-classical exact formula of functional increment. On the base of this formula the original problem is reduced to a problem of optimal control by a system of ordinary dierential equations.

Key words: optimization, optimal control, canonical hyperbolic equations.

ОПТИМАЛЬНЫЕ ФОРМЫ ТЕЛ, РАЗРУШАЮЩИХСЯ ЗА СЧЕТ РАДИАЦИОННОГО НАГРЕВА ПРИ ДВИЖЕНИИ В АТМОСФЕРЕ М.А. Аргучинцева, Э.В. Мошняков Иркутский государственный университет, Иркутск e-mail: marguch@math.isu.ru, savagem@mail.ru Аннотация. В данной статье исследована задача о нахождении оптимальной начальной формы аблирующего осесимметричного тела, движущегося по баллистической траектории и обладающего минимальным суммарным радиационным нагревом поверхности.

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

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

Нужно отметить, что унос массы теплозащитного покрытия приводит к изменению первоначальной формы тела, его аэродинамических коэффициентов и, в конечном счете, влияет на траекторию движения тела. Поэтому важно совместно решать задачи обтекания и нагрева тела с учетом абляции поверхности и движения тела по траектории. В наиболее общей постановке проблема поиска оптимальной начальной формы тела, обгарающего в атмосфере планеты, представляет собой сложную задачу оптимального управления систе мой нестационарных уравнений Навье-Стокса с управляемыми границами. Исследование этой задачи сопряжено со значительными трудностями. Во-первых, на каждой итерации оптимизационного метода необходимо решать сложную краевую задачу для системы нели нейных уравнений в частных производных. Во-вторых, пока еще не создан достаточно эффективный метод исследования нестационарных задач аэродинамики с учетом тепло обмена и абляции. В-третьих, поставленная задача является неклассической задачей опти мального управления с управляемыми границами, в качестве которых выступают форма тела и связанный с ней отход ударной волны. Таким образом, область изменения пере менных не задана, а определяется из решения задачи. В связи с вышеперечисленными Работа выполнена при финансовой поддержке РФФИ (проект N 05-01-00187) и программы "Универ ситеты России"(проект ур.03.01.064) трудностями проблема детального анализа общей задачи нахождения оптимальных форм тел вдоль траектории движения остается открытой.

В настоящее время в научной литературе рассматриваются различные упрощения ис ходной постановки задачи. Так, в работах [1,2] исследовалось уравнение абляции, когда унос массы происходил под воздействием конвективного нагрева. Соответствующие ре шения для стационарных форм аблирующих тел получены в работах [3-5]. Исследованию уноса массы осесимметричных и пространственных тел за счет радиационного теплового потока в атмосферах Земли и Юпитера при неизменных условиях в набегающем пото ке посвящены работы [6-9]. В [8] рассмотрена постановка сопряженной задачи, которая включает уравнения движения тела переменной формы и уравнения радиационной газо вой динамики для химически равновесных течений невязкого газа в ударном слое. По лученные в [8] результаты показали, что изменение формы аппарата в полете заметно влияет на траекторные параметры, тепловые потоки и характеристики теплозащиты. В работе [9] была решена задача об оптимальной стационарной форме аблирующего тела с точки зрения минимума уноса массы при неизменных условиях в набегающем потоке.

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

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

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

2) стационарную задачу обтекания тела в заданной точке траектории с учетом ра диационного теплообмена. Асимптотическое решение данной задачи позволило получить функционалы для радиационного теплового потока [10] и волнового сопротивления [11], явно зависящие от формы тела и параметров обтекания.

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

1. Постановка задачи Рассмотрим гиперзвуковое движение осесимметричного тела в атмосфере планеты вдоль баллистической траектории [8] dv = v 2 CB S + M g sin, M dt d g v = cos ( ), dt v RE dH = v sin. (1) dt Здесь t – время (t [0, T ]);

v = v(t), CB = CB (t), M = M (t), R = R(t), S = S(t) (S = R2 )– соответственно скорость, коэффициент волнового сопротивления, масса, радиус и площадь донного сечения тела;

= (t) – плотность газа на высоте H = H(t);

= (t) – угол наклона траектории движения;

g – ускорение свободного падения;

RE – средний радиус планеты. Для изотермической атмосферы имеем = 0 exp(H), где 0 – плотность атмосферы на уровне поверхности планеты;

1 – шкала высот для плотности.

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

yt qR y y =, yt =, yx =, (2) Hef T t x 1 + yx где x – координата вдоль оси тела;

y = y(x, t) – функция формы тела;

qR = qR (x, t) - величина локального радиационного теплового потока в заданную точку поверхности тела;

Hef - эффективная энтальпия разрушения теплозащитного покрытия;

T – плотность тела.

Уравнения (1), (2) решаются при начальных условиях H(0) = H0, v(0) = v0, (0) = 0 ;

(3) y(x, 0) = y0 (x), y0 (0) = 0, y0 (L) = R0, (4) где H0, v0, 0 – заданные высота, скорость и угол входа в атмосферу;

y0 (x), L, R0 – на чальная форма тела, ее длина и радиус донного сечения.

Выражения для коэффициента волнового сопротивления и локального радиационного теплового потока осесимметричного тела получены в работах [10,11] из асимптотического исследования задач стационарного гиперзвукового обтекания тел в заданной точке траек тории с учетом теплообмена L 4 y(x, t) yx (x, t) dx CB = 2, (5) R 1 + yx (x, t) x x v3 1 y(s, t) ys +11 (s, t) 2N qR = CR (1 + ys (s, t))(2N +9)/ 2 y(x, t) x N + N + x ys +8 (s, t) 2N 1 + b 1 + ys1 (s1, t) ds1 ds. (6) ys (s, t))(2N +7)/ (1 + s Здесь x0 = x0 (t) – координата критической точки тела (в дальнейшем будем предпола гать x0 (t) 0);

CR0 = CR0 (t) – коэффициент радиационного нагрева в критической точке тела;

N – параметр излучения газа [10]. Для вычисления коэффициента радиационного теплообмена в критической точке затупленного тела можно использовать известные эм пирические зависимости [10] для различных составов атмосфер планет. Например, для воздуха при 9 км/с v 20 км/с;

107 г/см3 104 г/см3 ;

0.3 м R 40 м имеем = 0.00344(0.1v)6 ( 106 )0.35 R0.5, CR0 =, b = 2(n + 4), 2 + 8.95/ Поставим следующую оптимизационную задачу: в классе гладких функций, удовле творяющих граничным условиям (4), найти начальную форму аблирующего тела y 0 (x) = y(x, 0), которая подчиняется системе уравнений (1)–(3) и минимизирует функционал сум марного радиационного нагрева T L QR = 2 y(x, t) (1 + yx (x, t))qR (x, t) dx dt. (8) 0 x Таким образом, поставленная задача является задачей оптимального управления на чальными условиями системы, состоящей из нелинейного интегро-дифференциального уравнения абляции (2) и системы обыкновенных дифференциальных уравнений движения тела по баллистической траектории (1). Для ее решения использовался численный метод локальных вариаций [12]. Опишем общую схему метода.

2. Метод решения Для дальнейших рассуждений функцию формы тела y(x, t) удобно переписать в без размерном виде (, t) ( [0, 1], t [0, T ]):

= y/R, = x/L, = R/L, = yx /, t = yt /R.

Тогда искомая начальная форма тела 0 () = y0 (x)/R будет удовлетворять граничным условиям 0 (0) = 0, 0 (1) = 1. На плоскости (, t) ведем равномерную сетку узлов ti = i t, t = T /n, i = 0, 1,..., n;

j = j, = 1/m, j = 0, 1,..., m, где соотношение между количеством узлов m и n выбиралось из условия выполнения соответствующего критерия Куранта для конечно-разностной схемы уравнения абляции (2). Первоначальное число узлов по временной переменной n = 80 затем уточнялось в процедуре автоматического выбора шага сетки t метода Рунге-Кутта интегрирования системы уравнений движения (1).

В качестве начального приближения искомого решения 0 () задавалась одна из сле дующих зависимостей:

0 = 1/4, 0 = 1/2, 0 = 3/4, 0 =.

0 0 0 Как показали расчеты, все эти начальные приближения в процессе варьирования сходи лись к одной и той же оптимизирующей кривой, однако, в случае конического начального приближения 0 = скорость сходимости метода была значительно быстрее. Для вы бранного начального приближения по квадратурным формулам Симпсона производился подсчет интегралов CB (5), M (7), qR (6), стоящих в правых частях уравнений движения (1) и абляции (2). Далее, на каждом временном слое проводилось совместное интегри рование указанной системы интегро-дифференциальных уравнений (1)–(2). Нужно отме тить, что система уравнений движения (1) решалась методом Рунге-Кутта 4-го порядка, а для решения уравнения абляции (2) была применена явно-неявная разностная схема. На рис. 1 представлены примеры расчетов начальной и конечных форм аблирующих тел при v0 = {11;

15;

19} км/c, H0 = 80 км, 0 = /3, 0 = 1.225 г/м3, 1 = 0.149 км, RE = км, = 0.5, T = 2200 кг/м3, Hef = 2 107 Дж/кг. Из расчетов следует, что конечная форма тела при ограничениях (4) стремится к конусу.

В конце процедуры вычислялось значение суммарного радиационного нагрева поверх ности вдоль траектории полета QR (8). Затем задавался шаг варьирования h по ординате из условия [12]: 0, h/( )2 0, необходимого для сходимости метода ло кальных вариаций. Процесс варьирования ординат 0 (j ) заключался в следующем: зна 0 чения 0 (j ) последовательно заменялись на (0 (j ) ± h). Для новых точек находились приближенные значения функционалов массы, сопротивления и нагрева (5)–(7), интегри ровалась система уравнений (1)–(2), и находилось значение целевого функционала QR (8).

Если значение функционала QR уменьшалось, то в таблицу решений записывалось соот ветствующее значение ординаты (0 (j ) ± h). Таким образом, проводился полный перебор 0 точек (j, 0 (j ) ± h), j = 0, 1,..., m, и находилось первое приближение экстремали {0 }, целевой функционал которого не больше вычисленного на предыдущем приближении.

В процессе счета итерации с фиксированными h и проводились до тех пор, k пока в результате очередной итерации не изменялось ни одно из значений {0 (j )}.

Затем шаг h уменьшался вдвое, и процесс варьирования продолжался. Уменьшение h при фиксированном проводилось до тех пор пока не удовлетворялось неравенство h 1010. Затем рассматривалось решение задачи на более мелкой сетке, т.е. число разбиений m удваивалось. Функция = () во вновь образованных точках находилось интерполяцией по соседним узлам. Процесс варьирования продолжался до тех пор, пока не выполнялось неравенство 105.

3. Анализ результатов расчетов Была проведена серия расчетов оптимальных начальных форм осесимметричных тел в широком спектре условий входа в атмосферу Земли и заданных геометрических ха рактеристиках тела. Например, на рис.2 представлены оптимальные начальные формы аблирующих тел при следующих параметрах входа в атмосферу Земли: v0 = 15 км/c, H0 = 80 км, 0 = /3, 0 = 1.225 г/м3, 1 = 0.149 км, RE = 6371 км. Соответствующие заданные геометрические и теплофизические характеристики тел: = {1. 0.5;

2. 0.8}, T = 2200 кг/м3, Hef = 2 107 Дж/кг.

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

Рис. Рис. Список литературы [1] В.В. Лунев Некоторые свойства уравнения абляции. Изв. АН СССР. МЖГ, 1977, № 3, с.96-102.

[2] В.В. Знаменский Численное решение уравнения уноса. Изв. АН СССР. МЖГ, 1978, № 2, с.147-154.

[3] В.Г. Воронкин, В.В. Лунев, А.Н. Никулин О стационарной форме тел при их разру шении за счет аэродинамического нагрева. Изв. АН СССР. МЖГ, 1978, № 2, с.138-146.

[4] И.Н. Мурзинов Аналитическое решение о стационарных формах тел в условиях аб ляции. В кн.: Газовая и волновая динамика. М.: Изд-во Моск. ун-та, 1979.

[5] В.Г. Коняев Аналитическое исследование изменения формы аблирующих тел при их движении в атмосфере со сверхзвуковыми скоростями. Уч. записки ЦАГИ, 1974, № 6.



Pages:     | 1 || 3 | 4 |   ...   | 6 |
 





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

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