Решение задачи оптимизации многоступенчатых ракет

cl-jynieyqxwxyasvn9tfcst0jg.png

Введение

Методы нелинейной оптимизации широко применяются при проектировании машин и механизмов. Указанные методы применяются и в ракетостроении, например, для оптимизации многоступенчатых ракет [1].

Многоступенчатая ракета — это аппарат, в котором части конструкции отделяются во время полета, придавая оставшейся части ракеты дополнительную скорость. Трёхступенчатая ракета схематически показана на рисунке.

lwfpkcmed2d0lhdigvold_76qgy.png

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

Мы рассмотрим две задачи в предположении, что коэффициент kaebuvntumpsdwf4-ff8ivuiubu.png и скорость реактивной струи Cn постоянны на каждой ступени, однако на разных ступенях могут принимать различные значения. В обеих задачах в качестве целевой функции принят коэффициент полезной нагрузки ракеты G, который необходимо минимизировать.

Характеристики многоступенчатой ракеты можно описать двумя уравнениями. Первое уравнение для коэффициента полезной нагрузки ракеты:

r2lvchdx5w7_ngo-mhjufelbg7e.png

где: W1– полезный вес ракеты ; WN –начальный вес ракеты до отделения ступеней.

Второе уравнение для характеристики полезной нагрузки многоступенчатой ракеты, которое можно записать как выражение для еѐ идеальной конечной скорости (формула К.Э. Циолковского):

oue9byv7kef7v0swkxve9objmmk.png

где:

zruauxw-fvhl4dqt4oeuefeskpg.png

Рассмотрим 4-ступенчатую ракету, движущуюся в отсутствии гравитационного поля с заданными значениями идеальной конечной скорости Vc, начального веса ракеты W4, скорости реактивной струи Cn и коэффициентов Ɛn=Ɛ.

Постановка задач

Первая задача — найти такое распределение веса по ступеням W3, W2, чтобы минимизировать коэффициент полезной нагрузки ракеты G=W4/W1 при заданной скорости Vc.

Вторая задача — найти такое распределение веса по ступеням W3, W2, чтобы максимизировать идеальную конечную скорость Vc при заданном коэффициенте полезной нагрузки ракеты G .

Рассмотрим первую задачу минимизации коэффициента полезной нагрузки при заданной скорости ракеты Vc =1400 м/c и скорости реактивной струи Cn = C =2000 м/с.

В качестве управляющих переменных выбираем распределение веса по ступеням: W3, W2 при заданных значениях начального веса ракеты W4 и коэффициента влияния нагрузки Ɛn=Ɛ=0.3

Представим формулу для вычисления заданной скорости в виде:

6xu4c2isnnhz3ztbcf6tsih8-rk.png

где

gbxabc9nu2gnf18vfdlrfdobcdq.png

Отсюда следует, что:

7utmzgdwxs2tjrb0cl5djxjh7ue.png

Решение первой задачи на Python

Python является свободно распространяемым языком программирования высокого уровня. До настоящего времени задачи оптимизации решались в лицензионных математических пакетах Maple, Mathcad, Matlab и других.

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

Привожу листинг первой части программы, включая функцию цели:

Определение функции цели
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import time
C=2000 # скорость реактивной струи
e=0.3# коэффициент нагрузки
VCzad=1400# идеальная конечная скорость ракеты
w4=4000# начальная масса ракеты
def f(a):# вспомогательная функция
                return C*np.log(a/(1+e*a))
def VC2(w3,w2):#скорость после отделения первой ступени
                return VCzad-f(w4/w3)-f(w3/w2)
def w1(w3,w2):# функция для расчёта массы первой ступени
                return (w2/(np.e**(VC2(w3,w2)/C)))-e*w2
start = time.time()# начало отсчета времени работы блока оптимизации
def fun1(x,y):#функция цели с новыми переменными (x=w3, y=w2) в форме удобной для расчетов
         return 4000/(-0.3*y+y*np.e**(np.log(x/(0.3*x+y))+np.log(4000/(1200+x))-0.7))

Построим контурный график функции цели.

Листинг графика функции цели
#!/usr/bin/env python
#coding=utf8
import numpy as np 
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D 
def fun1(x,y):#функция цели 
         return 4000/(-0.3*y+y*np.e**(np.log(x/(0.3*x+y))+np.log(4000/(1200+x))-0.7))
x= np.linspace(500,5000,10000)
y = np.linspace(500,5000,10000)
z=fun1(x,y)
fig = plt.figure('Контурный график функции цели')
ax= Axes3D(fig)
ax.plot(x,y,z,linewidth=1)
plt.show()

Получим:

edzi0wndjxlbash8j0oqw9po6qy.png

Сложность функции цели состоит в резких изменениях в районе приемлемых значений переменных, поэтому исследуем библиотеку from scipy. optimize import minimize.

Обычно функцию цели записывают в виде [2]: fun = lambda x: (x[0] — 1)**2 + (x[1] — 2.5)**2 (пример). Использование такой записи в решаемой задачи усложнит анализ целевой функции. Для упрощения формы записи целевой функции, применим вспомогательную функцию fun2(x):

def fun2(x):
   return fun1(*x)

Используем пример из документации [2] для решателя SLSQP:

import numpy as np
from scipy.optimize import minimize
import time
start = time.time()
def fun1(x,y):
         return 4000/(-0.3*y+y*np.e**(np.log(x/(0.3*x+y))+np.log(4000/(1200+x))-0.7))
def fun2(x):
   return fun1(*x)
x0 =[3000,2000] # начальная точка поиска решения
res = minimize(fun2, x0, method='SLSQP', bounds=None)
stop = time.time()
print ("Время оптимизации :",round(stop-start,3))
print(res)

Получаем следующий результат:

Время оптимизации: 0.0
fun: 8.8345830746128
jac: array ([ 7.71641731e-04, -4.12464142e-05])
message: 'Optimization terminated successfully.'
nfev: 48
nit: 12
njev: 12
status: 0
success: True
x: array ([ 2953.94929027, 1149.31138351])

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

import numpy as np
from scipy.optimize import minimize
import time
start = time.time()
def fun1(x,y):
         return 4000/(-0.3*y+y*np.e**(np.log(x/(0.3*x+y))+np.log(4000/(1200+x))-0.7))
def fun2(x):
   return fun1(*x)
x0 =[3000,2000] # начальная точка поиска решения
res = minimize(fun2, x0))
stop = time.time()
print ("Время оптимизации :",round(stop-start,3))
print(res)

Получаем следующий результат:

Время оптимизации: 0.0
fun: 8.402280475718243
hess_inv: array ([[ 911166.17713172, 192191.3238931 ],
[ 192191.3238931, 194600.77394409]])
jac: array ([ -1.19209290e-07, 1.31130219e-06])
message: 'Optimization terminated successfully.'
nfev: 152
nit: 28
njev: 38
status: 0
success: True
x: array ([ 1967.29278561, 967.91321904])

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

Листинг решения первой задачи оптимизации
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import time
C=2000 # скорость реактивной струи
e=0.3# коэффициент нагрузки
VCzad=1400# идеальная конечная скорость ракеты
w4=4000# начальная масса ракеты
def f(a):# вспомогательная функция
                return C*np.log(a/(1+e*a))
def VC2(w3,w2):#скорость после отделения первой ступени
                return VCzad-f(w4/w3)-f(w3/w2)
def w1(w3,w2):# функция для расчёта массы первой ступени
                return (w2/(np.e**(VC2(w3,w2)/C)))-e*w2
start = time.time()# начало отсчёта времени работы блока оптимизации
def fun1(x,y):#функция цели с новыми переменными (x=w3, y=w2) в форме удобной для расчётов
         return 4000/(-0.3*y+y*np.e**(np.log(x/(0.3*x+y))+np.log(4000/(1200+x))-0.7))
def fun2(x):
   return fun1(*x)
x0 =[3000,2000]
res = minimize(fun2, x0)
stop = time.time()
print ("Время работы оптимизатора :",round(stop-start,3))
print ("Расчётное значение функции цели :",round(res['fun'],3))
w2=round(res['x'][1],3)
w3=round(res['x'][0],3)
w1=round(w1(w3,w2),3)
v1=round(f(w2/w1),3)
v2=round(f(w2/w1)+f(w3/w2),3)
v3=round(f(w2/w1)+f(w3/w2)+f(w4/w3),3)
print("Изменение веса ракеты -:%s;%s;%s;%s "%(w4,w3,w2,w1))
print("Изменение скорости ракеты -:%s;%s;%s;%s "%(0,v1,v2,v3))
VES=[w4,w3,w2,w1]
VC=[0,v1,v2,v3]
plt.title("Распределение веса ракеты при минимальной нагрузке \n для заданной скорости")
plt.plot(VES,VC,'o')
plt.plot(VES,VC,'r')
plt.xlabel("Вес ракеты")
plt.ylabel("Скорость ракеты")
plt.grid(True)          
plt.show()

Получаем следующий результат:

Время работы оптимизатора: 0.0
Расчётное значение функции цели: 8.402
Изменение веса ракеты -:4000;1967.293;967.913;476.061
Изменение скорости ракеты -:0;466.785;933.166;1400.001

lydokm4hs7wpnocgy0dvvddff78.png

Вес ракеты распределен так (4000.0; 1967.0; 967.9; 476.1), что соответствует нарастанию скорости равными долями (после отделения каждой ступени) в среднем на 460 м/с.

Используем ту же формулу для конечной скорости ракеты, что и в предыдущей задаче:

j77pejs2tpovsrwwlj_xli8flq8.png

Так как W4 и W1 заданы, то распределением веса ракеты по ступеням W3 и W2 можно максимизировать идеальную конечную скорость Vc при заданном коэффициенте полезной нагрузки ракеты G.

Решение второй задачи на Python

С учётом ранее изложенного, запишем часть листинга до функции цели:

Листинг для определения функции цели
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import time
C=2000 # скорость реактивной струи
e=0.4# коэффициент влияния изменения нагрузки
w4zad=4000# масса ракеты
w1zad=200# масса первой ступени
def f(a):
         return np.log(a/(1+e*a))
def Vc1(w3):
         return C*f(w4zad/w3)
def Vc2(w3,w2):
         return C*f(w3/w2)
def Vc3(w2):
         return C*f(w2/w1zad)
start = time.time()
def fun1(x, y):
         return  -(C*np.log((w4zad/y)/(1+e*(w4zad/y)))+ C*np.log((y/x)/(1+e*(y/x)))+C* np.log((x/w1zad)/(1+e*(x/w1zad))))

Построим контурный график функции цели:

Листинг контурного графика
import numpy as np 
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
C=2000 #скорость реактивной струи
e=0.4#коэффициент влияния изменения нагрузки
w4zad=4000# масса ракеты
w1zad=200# масса первой ступени
def fun1(x, y):
         return  -(C*np.log((w4zad/y)/(1+e*(w4zad/y)))+ C*np.log((y/x)/(1+e*(y/x)))+C* np.log((x/w1zad)/(1+e*(x/w1zad))))
x= np.linspace(500,5000,10000)
y = np.linspace(500,5000,10000)
z=fun1(x,y)
fig = plt.figure('Контурный график функции цели')
ax= Axes3D(fig)
ax.plot(x,y,z,linewidth=1)
plt.show()  

Получим контурный график функции цели с отрицательным знаком для смены минимума на максимум:

brbh2pxptg3edycyn_h7c13j3da.png

Получили локальный максимум (минимум на графике), поэтому окончательный листинг представим в виде:

Листинг решения второй задачи оптимизации
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import time
C=2000 #скорость реактивной струи
e=0.4#коэффициент влияния изменения нагрузки
w4zad=4000# масса ракеты
w1zad=200# масса первой ступени
def f(a):
         return np.log(a/(1+e*a))
def Vc1(w3):
         return C*f(w4zad/w3)
def Vc2(w3,w2):
         return C*f(w3/w2)
def Vc3(w2):
         return C*f(w2/w1zad)
start = time.time()
def fun1(x, y):
         return  -(C*np.log((w4zad/y)/(1+e*(w4zad/y)))+ C*np.log((y/x)/(1+e*(y/x)))+C* np.log((x/w1zad)/(1+e*(x/w1zad))))
def fun2(x):
   return fun1(*x)
x0 = (1000,2000)
res = minimize(fun2, x0)
stop = time.time()
print ("Время работы оптимизатора:",round(stop-start,3))
print ("Расчётное значение функции цели :",round(res['fun'],3))
w2=round(res['x'][0],3)
w3=round(res['x'][1],3)
v1=round(Vc1(w3),3)
v2=round(Vc1(w3)+Vc2(w3,w2),3)
v3=round(Vc1(w3)+Vc2(w3,w2)+Vc3(w2),2)
print("Изменение веса ракеты -:%s;%s;%s;%s "%(w4zad,w3,w2,w1zad))
print("Изменение скорости ракеты -:%s;%s;%s;%s "%(0,v1,v2,v3))
VES=[w4zad,w3,w2,w1zad]
VC=[0,v1,v2,v3]
plt.title("Распределение веса ракеты при максимальной скорости\n для заданной нагрузки")
plt.plot(VES,VC,'o')
plt.plot(VES,VC,'r')
plt.xlabel("Вес ракеты")
plt.ylabel("Скорость ракеты")
plt.grid(True)          
plt.show()

В результате получим:

Время работы оптимизатора: 0.016
Расчётное значение функции цели: -1580.644
Изменение веса ракеты -:4000;1473.625;542.889;200
Изменение скорости ракеты -:0;526.873;1053.753;1580.64

tlkl2kxzpai1tibd24awmoywu0c.png

Веса ракеты распределяются так:4000;1473.625;542.889;200, что соответствует нарастанию скорости равными долями (после отделения каждой ступени) на 526.9 м/с.

Выводы

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

Всем спасибо за внимание! И главное, чтобы многоступенчатая ракета с не полезной для нас нагрузкой НИКОГДА к нам не прилетела.

Ссылки

  1. Решение задач нелинейного программирования
  2. SciPy. Optimize.minimize

© Habrahabr.ru