Применение метода множителей для решения задачи баллистического проектирования ракеты-носителя Текст научной статьи по специальности «Математика»

УДК 629.7.01

Применение метода множителей для решения задачи баллистического проектирования

ракеты-носителя

© К.П. Баслык, Н.Н. Генералов, Б.Г. Кулешов МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

Рассмотрена проблема баллистического проектирования трехступенчатой ракеты-носителя. Задача нелинейного программирования сформулирована как определение минимума функции цели при наличии ограничений. В качестве целевой функции выбрана величина, обратная относительной массе полезного груза, выводимого на заданную орбиту. Ограничения заданы в виде параметров конца активного участка траектории ракеты-носителя. При решении задачи на условный экстремум применен метод множителей в сочетании с методами покоординатного спуска и золотого сечения для решения задач безусловной оптимизации в одно- и многомерном случаях. Для численного интегрирования уравнений движения применяется метод Рунге — Кутты 4-го порядка точности. Предложен алгоритм проектировочного расчета, разработана оригинальная программа на языке Паскаль, проведен ряд тестовых расчетов.

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

Баллистический анализ трехступенчатой ракеты-носителя (РН) в предлагаемой в статье формулировке представляет собой решение задачи нелинейного программирования, а именно, нахождение минимума целевой функции, обратной относительной массе полезного груза (ПГ) ракеты:

= (1) - ПГ ГТ-

1=1

при заданных ограничениях вида

V = УА, Нк = ИА; = 0. (2)

где ¥к, Ик, 3 к — скорость и высота полета, а также угол наклона к местному горизонту в конце активного участка траектории (АУТ), км/с, км, град соответственно.

Значение целевой функции 1/ - ПГ (см. формулу (1)) однозначно определяется неизвестными относительными конечными массами ступеней ракеты -кг, а также массовыми коэффициентами отсеков и

агрегатов ракетных блоков, значения которых считаются известными. Помимо этого по ступеням ракеты задаются энергетические характеристики маршевых двигательных установок (ДУ) ступеней — удельные импульсы тяги и начальные нагрузки на тягу. Относительная масса ПГ 1-й ступени ракеты [1]

М ПГ г = Д'М* к г _ Вг, (3)

где Лг, Вг — коэффициенты, определяемые соответственно по формулам

Л = 1 + а£ г ; в = М£ г +УДУ г/ У0 г + а£ г . (4)

1 + МПЕРг 1 + МПЕР г 1 + МПЕР г

У ду < = , (5)

Р0 г

МдУ г — масса конструкции ДУ, кг; g0 = 9,81-10-3 км/с2 — ускорение свободного падения на Земле; Р0 г — начальная тяга ДУ ступени, кН; у0 — стартовая нагрузка на тягу (для второй и третьей ступени используют нагрузку уп на тягу в пустоте).

Относительные массы (коэффициенты а и м) отсеков и агрегатов вычисляют по формулам

М£ г = М ПО г +М ПУ г +М ХО г +М ОУ г. (6)

а£г = аТО г + аСП г ат.ост г. (7)

Здесь

Мг (8)

М' = М^; (8>

а. = Мт-; (9)

М т г

Мг — масса агрегата или отсека г-й ступени; М0г — начальная масса г-й ступени; Мтг — начальная масса топлива, т. е. расчетной характеристикой для м является начальная масса М0 г г-й ступени, а для а — запас топлива М т г.

В формулах (4)-(9) индексы обозначают: ПО — приборный отсек; ПУ — приборы системы управления; ХО — хвостовой отсек; ОУ — органы управления; ПЕР — переходной отсек; ТО — топливный отсек; СП — система подачи; т.ост — топливные остатки.

Параметры конца АУТ вычисляют в результате интегрирования уравнений движения, записанных в виде [1-3]

ёУ ё |

1 уд 0 а э ^ р н

1уд 0 ( а э 1)

I

Ро

соб а +

СхЧУ 01 уд о ^ Я

IР ёН

Я о

V о I уд оБш 3;

V

=--V о I уд оэт 3;

ё| Яо

ё 1 ё |

V

(Я® + Н )

-V о I уд о соб 3;

(10)

(11) (12)

ё3 ё |

У

1 ул оа э Рн 1 уд о (аэ - 1) ^

+ о1 уд о +

I Р м

1 уд 0а э

I

Я

Ро

I

Бт а +

V2

Я о (®+ Н )я о

V оIуд о соб 3!

(13)

Здесь (индекс г в дальнейшем не указываем) V, Н — скорость, км/с, и высота полета, км, соответственно; 7*уд о — удельный импульс тяги, км/с; аэ — коэффициент увеличения удельного импульса в пустоте; РН, Ро — давление окружающей среды на высоте полета и на Земле; v0 — стартовая нагрузка на тягу (для второй и третьей ступени используют vп — нагрузку на тягу в пустоте); Сх, Су — аэродинамические коэффициенты лобового сопротивления и подъемной силы соответственно; я0, я — ускорение свободного падения на Земле и на высоте полета соответственно, км/с2; а — угол атаки; q = pV 2/2 — скоростной напор, значение которого следует вычислять в таких единицах, чтобы в результате второе слагаемое в формуле (13) измери-лось бы в км/с; Рм — стартовая нагрузка на мидель, принятая равной 0,1 МПа; 1 — полярный угол (координата); Я® — радиус Земли, км.

Аэродинамические коэффициенты лобового сопротивления и подъемной силы с учетом малости углов атаки на атмосферном участке полета могут быть представлены как

Сх (М, а) = Сх 1 соб а + Су 1 Бт а =

Су (М, а) = — Сх1 Бт а + Су 1 соб а = (Суа - С°1)

бш а^а; соб а «1;

Су 1 = СУ1а;

С = С о

^Х 1 — ^Х 1

= С 0 + С аа 2 _ С0 •

а,

(14)

где

С01 = Сх! (М,0);

С а = 1 ^-У1 _ '

5а

(15)

При интегрировании уравнений движения (8)-(10) коэффициенты в формулах (15) принимают следующими [4]:

С о =

^ х 1 —

0,29 М - 0,51 0,5

при 0 < М < 0,8; при 0,8 < М < 1,068;

(16)

0,091 + -

М

при

М > 1,068;

С а _ ^ У1 —

2,8 при 0 < М < 0,25; 2,8 + 0,447(М-0,25) при 0,25 < М < 1,1; 3,18-0,660(М-1,1) при 1,1 < М < 1,6; 2,85 + 0,350(М-1,6) при 1,6 < М < 3,6; 3,55 при М > 3,6.

(17)

Для моделирования выведения ПГ была использована двухпара-метрическая программа тангажа РН, формируемая из четырех последовательно отрабатываемых участков [3, 5]:

вертикальный полет РН до м = 0,96 на участке первой ступени; участок программного разворота по углу атаки на первой ступени; участок гравитационного разворота до конца работы первой ступени;

участок программного разворота по углу тангажа с начала работы второй ступени и до конца АУТ.

Программные законы изменения угла атаки на первой ступени и угла тангажа на второй и третьей ступенях задают в виде [3]

апр ^) = -4атеа(-) (1 - еа(-));

ф пр ( ( ) = ф 0 +ф пр ( ( - (2 ) ,

(18) (19)

где а т — максимальный угол атаки; — время окончания участка вертикального полета; ф0, Ь — угол тангажа и время окончания работы первой ступени соответственно.

Таким образом, в качестве параметров программы полета выступают величины а т и фпр. Время полета (в пределах работы ступени) в формулах (18), (19) определяется выражением

а=0

t = (1 -^)vo (20)

g o

Значение коэффициента a в формуле (18) при расчетах принимают равным 0,30...0,40 с-1. Это обеспечивает программный разворот ракеты на первой ступени продолжительностью 20.30 с, а также пренебрежимо малые углы атаки при достижении ракетой скорости порядка 0,250 км/с.

Остановимся на методах решения задачи условной оптимизации. Для задачи с ограничениями типа равенств (2) математическая формулировка имеет вид [6-11]:

X*: f (X*min; h (X * ) = 0,

где X* — вектор неизвестных; f (X*) — целевая функция, определяемая формулой (1); h(X ) — вектор ограничений (см. равенства (2)).

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

F1 (X, X) = f (X ) + X?h (X ), (21)

где X — вектор множителей Лагранжа, т. е. определяют вектор методами непрямого поиска [12, 13], например методом Ньютона:

(X*T X*T )T : VF (X*, X*) = 0. (22)

Один из недостатков этого подхода состоит в необходимости вычисления матрицы Гессе, содержащей вторые производные от функции F1 (см. формулу (21)), которые, впрочем как и первые, могут быть определены только численно (см. выражения (10)-(13)). Другое неудобство использования этого метода связано с тем, что для обеспечения устойчивости вычислений по схеме Ньютона и достаточно быстрой сходимости метода необходимо задать хорошее начальное

приближение к вектору (X*T X*T )T, что практически невозможно

из-за неизвестного заранее значения вектора X.

Применение метода штрафных функций [6-11, 13], в частности простого квадратичного штрафа, состоит в том, что строят последовательность параметров штрафа {ck j, ck ^да и для каждой итерации (фиксированного значения cN) решают задачу безусловной минимизации вида

X" : F2 (X",^) = f (XN) + 0,5cN |и(XN)) ^ шт.

2

Критерием окончания итерационного процесса является достижение

где в — заданная точность решения.

Доказательство сходимости метода штрафа приведено, например, в работах [7, 8]. Однако в них отмечено, что при больших значениях параметра штрафа вспомогательная задача безусловной оптимизации становится плохо обусловленной. В частности, для ее решения гарантированно неприменимы методы наискорейшего спуска.

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

где значения параметра штрафа c и вектора множителей Лагранжа X пересчитываются на каждой итерации.

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

X N : F2 (X N, XN, ^ ) = f (XN ) + XNТИ (XN ) + 0,5^ |и (XN )|2 ^ шт,

2

F3 (X, X ь X 2, X з, е ) = —

П ПГ г

1=1

относительно вектора неизвестных

Х — [а т , Д кЪ Д к 25 Д к35 ф пр ] •

(24)

Здесь Х1з Х2, Х3 — неопределенные множители Лагранжа.

Алгоритм решения задачи следующий [6]:

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

2) задание начальных значений множителей Лагранжа X0, г = 1, 2, 3, а также параметра штрафа с ;

3) решение задачи многомерной оптимизации для фиксированных значений множителей Лагранжа и параметра штрафа с заданной точностью, т. е. определения минимума функции ^з(Х, Х1, Х2, Х3, с) (см. формулы (23), (24));

4) вычисление значений множителей Лагранжа:

X N+1 — X N + с'

X N+1 — х N + с>

+ сн ( V?

V Vа

+ сы ( Н м ±± к

V Н а

X N + ^ 3 N

-1

-1

5) определение параметра штрафа по формуле

см+1 — р^;

(25)

6) переход к п. 3, если точность по параметрам Ук, Нк и 3 к полученного решения задачи недостаточная.

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

воспользуемся формулой Циолковского. Тогда

( \

д к 1 — д к — ехр

1,22Уа

3

^ \а э г1 уд 0г

V 1—1

г = 1, 2, 3.

Начальные значения для параметров программы полета (см. формулы (18), (19)) принимают равными ат « 0,5° и фпр « 0,005 °/с. Та-

кие небольшие значения исключают излишнее заваливание ракеты на начальном участке траектории.

Численные эксперименты показали, что для широкого диапазона значений исходных данных (значений у0г- и /уд0г) при решении поставленной задачи (см. формулы (23), (24)) алгоритм метода множителей устойчиво работает для следующих начальных значений множителей Лагранжа, параметра штрафа и коэффициента Р из формулы (25):

X 0 =Х2 =Х0 =0;

е0 = 800; (26)

Р = 10.

Задачу многомерной оптимизации при фиксированном значении множителей Лагранжа и параметра штрафа (см. п. 3 алгоритма) предлагается решать методом покоординатного спуска [12] с коррекцией на каждой итерации интервала поиска. Этот метод безусловной оптимизации позволяет гарантированно избежать выхода за границы интервала поиска решения, а также не требует вычисления производных от модифицированной функции Лагранжа. Для каждой компоненты вектора X (см. формулу (24)) задача одномерной оптимизации решается методом золотого сечения. После выполнения нескольких итераций (в реализованной программе число таких итераций равно 5) для уменьшения вычислительных затрат на каждом шаге проводят коррекцию интервала поиска решения. Границы нового интервала [а?, bN] вычисляют по формулам

а? = х? 1 -5?-1;

Ь? = х? 1 + 5?-1,

где 5? 1 = |х? 1 - х?-2|, х? — ¡-я компонента из вектора X, определенная на ?-й итерации.

С увеличением параметра штрафа е может наблюдаться рост числа обусловленности задачи одномерной оптимизации (см. [12]). Для снижения влияния этого негативного фактора применен следующий искусственный прием. Требования по точности получаемых решений ужесточают последовательно с увеличением числа итераций ? (рис. 1). При этом для каждого нового порогового значения точности решения Еps_i обнуляются начальные значения множителей Лагранжа и устанавливается начальное значение параметра штрафа (см. формулу (26)).

О

Еря_ 1

Еря 2 Д |

у у V V V й

Рис. 1. Требования к точности решения задачи

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

б, =

3 Л, N

Ц" . _Ц N-1 ^ 2 ^ _а ^-1 ^ 2 Гф - (0 ^ 2

г^к г г^к г , т т , т^пр т^пр

Ц

N-1

10 V180

5я/180

и точность решения по параметрам конца АУТ

Б 2 =

г и

V VА

_ 11 +

Г нм

11 к

V Н А

-1| + )2

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

Первый расчет представляет собой баллистическое проектирование РН, близкой по своим характеристикам к ракете «Протон». Масса выводимого ПГ на круговую орбиту высотой 250 км составляет 23 т. Точности решения задачи приведены в табл. 1, а исходные данные и результаты расчета — в табл. 2.

Таблица 1

Требуемая точность получаемых решений для проектируемой РН

г=1

№ По проектным параметрам РН По параметрам конца АУТ

1 0,1000 0,1000

2 0,0357 0,0455

3 0,0128 0,0207

4 0,0046 0,0094

5 0,0016 0,0043

Таблица 2

Исходные данные и результаты баллистического проектирования РН

Сту- Исходные данные Результаты расчета

пень 7уд0, км/с аэ Це М-пер Уду Га, км/с км а, с 1 Цк М0,т х, с км/с Як, км а ° й?фпр/й?/, °/с

I 2,845 1,16 0,65 0,015 0,040 0,0551 0,0130 7,759 250,0 0,300 0,3292 804,2 126,4 7,735 249,813 0,113 2,114 -0,121

II 3,159 1,0 0,90 0,016 0,030 0,0551 0,0112 7,759 250,0 0,300 0,4296 199,0 165,3 7,735 249,813 0,113 2,114 -0,121

III 3,159 1,0 1,00 0,018 0,0 0,0551 0,0123 7,759 250,0 0,300 0,3861 71,4 197,7 7,735 249,813 0,113 2,114 -0,121

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

3

^хар = —^ аэг1уд0г 1п М* к г = г =1

= -(1,16 • 2,8451п 0,3292 + 3,1591п0,4296 + 3,1591п 0,3861) =

= 9,342 км/с.

Затем составим функцию Лагранжа вида

^ (м к1, Мк2 , М к3, =

= ~3 + Х [ ^хар + (э11уд 0 1 1п М к 1 + аэ21 уд0 2 1п М к 2 + аэ31уд 0 3 1п М к3 ) ,

Пм ПГ г г =1

для которой методом Ньютона определяем точку экстремума.

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

относительные конечные массы ступеней мк1 = 0,3528; мк2 = = 0,3885; м к3= 0,3971;

начальные массы ступеней М0 1 = 797,505 т; М0 2 = 217,561 т; М0 3 = 68,939 т.

На рис. 2 и 3 представлены результаты расчета параметров траектории проектируемой РН.

Расчеты были также проведены для трехступенчатой РН, выводящей ПГ массой 3,0 т в точку конца АУТ со второй космической скоростью 11,3 км/с [1]. Для заданных исходных данных получены проектные параметры РН, выводящих ПГ на круговые орбиты высотой 250, 300 и 350 км.

Задачу решали для значений точности, приведенных в табл. 3. Исходные данные и результаты расчета для этого случая представлены в табл. 4

Анализ данных табл. 4 показывает хорошее соответствие результатов, полученных в расчетной программе и приведенных в [1].

V, км/с

100 200 300 400 т, с 0 в

100 200 300 400 т, с

Рис. 2. Зависимости скорости РН и коэффициента продольной перегрузки (а), угла наклона траектории и высоты полета (б), а также угла тангажа и угла атаки (в) от времени полета на АУТ

Eps_2

Eps_\

Eps г

0,12 0,08 0,04

0

Eps_\

0,002 ----

0 80 100 N

40

80

N

Рис. 3. Зависимость точности вычисления проектных параметров РН (Eps_1) и параметров конца АУТ (Eps_2) от числа итераций

Таблица 3

Требуемая точность получаемых решений для тестовых расчетов

№ Требуемая точность получаемых решений

По проектным параметрам РН По параметрам конца АУТ

1 0,1000 0,1000

2 0,0417 0,0455

3 0,0174 0,0207

4 0,0072 0,0094

5 0,0030 0,0043

Таблица 4

Исходные данные и результаты баллистического проектирования для тестовых расчетов

Ступень Исходные данные Результаты расчета при ЯА, км Данные работы [1]

/уд0 км/с аэ ЦЕ М-ПЕР ЯЕ Уду УА, км/с а, с 1 250 300 350 Цк М0,т

Цк М0,т Цк М0,т Цк М0,т

I 2,747 1,16 0,55 0,015 0,040 0,0551 0,0130 11,300 0,400 0,2337 434,4 0,2311 465,5 0,2409 503,3 0,2960 486,0

II 3,384 1,0 0,48 0,016 0,030 0,0551 0,0112 11,300 0,400 0,3145 63,8 0,3090 67,2 0,3096 77,7 0,2650 102,5

III 3,434 1,0 0,54 0,018 0,0 0,0551 0,0123 11,300 0,400 0,2840 14,7 0,2790 15,1 0,2532 18,5 0,2480 18,5

ЛИТЕРАТУРА

[1] Разумеев В.Ф. Выбор проектно-баллистических параметров ракет. Москва, Изд-во МВТУ, 1973, 154 с.

[2] Мишин В.П., Безвербый В.К., Панкратов Б.М., Щеверов Д.Н. Основы проектирования летательных аппаратов (транспортные системы). Мишин В.П., ред. Москва, Машиностроение, 1985, 360 с.

[3] Феодосьев В.И. Основы техники ракетного полета. Москва, Наука, 1979, 496 с.

[4] Варфоломеев В.И., Копытов М.И. Проектирование и испытания баллистических ракет. Москва, Воениздат, 1970, 392 с.

[5] Сердюк В.К. Проектирование средств выведения космических аппаратов. Москва, Машиностроение, 2009, 504 с.

[6] Бертсекас Д. Условная оптимизация и методы множителей Лагранжа. Москва, Радио и связь, 1987, 400 с.

[7] Зангвилл У. Нелинейное программирование. Единый подход. Москва, Советское радио, 1973, 312 с.

[8] Полак Э. Численные методы оптимизации. Единый подход. Москва, Мир, 1974, 376 с.

[9] Поляк Б.Т. Введение в оптимизацию. Москва, Наука, 1983, 384 с.

[10] Реклейтис Г., Рейвиндран А., Рэгсдел К. Оптимизация в технике. Т. 1. Москва, Мир, 1986, 350 с.

[11] Химмельблау Д. Прикладное нелинейное программирование. Москва, Мир, 1975, 536 с.

[12] Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров. Москва, Высшая школа, 1994, 544 с.

[13] Турчак Л.И. Основы численных методов. Москва, Наука, 1987, 318 с.

Статья поступила в редакцию 06.07.2013

Ссылку на эту статью просим оформлять следующим образом: Баслык К.П., Генералов Н.Н., Кулешов Б.Г. Применение метода множителей для решения задачи баллистического проектирования ракеты-носителя. Инженерный журнал: наука и инновации, 2013, вып. 7. URL: http://engjournal.ru/catalog/machin/rocket/849 .html

Баслык Константин Петрович окончил МГТУ им. Н.Э.Баумана в 1994 г. Канд. техн. наук, доцент кафедры «Космические аппараты и ракеты-носители» МГТУ им. Н. Э. Баумана. Автор 15 опубликованных научных трудов в области механики и численных методов расчета конструкций из композиционных материалов. e-mail: [email protected]

Генералов Николай Николаевич окончил МВТУ им. Н.Э. Баумана в 1967 г. Старший преподаватель кафедры «Космические аппараты и ракеты-носители» МГТУ им. Н.Э. Баумана. Автор 27 научных трудов в области механики и проектирования летательных аппаратов. e-mail: [email protected]

Кулешов Борис Георгиевич окончил МВТУ им. Н.Э. Баумана в 1960 г. Старший преподаватель кафедры «Космические аппараты и ракеты-носители» МГТУ им. Н.Э. Баумана. Автор 28 научных трудов в области механики и проектирования летательных аппаратов. e-mail: [email protected]