![](https://habrastorage.org/webt/cl/-j/yn/cl-jynieyqxwxyasvn9tfcst0jg.png)
Введение
Методы нелинейной оптимизации широко применяются при проектировании машин и механизмов. Указанные методы применяются и в ракетостроении, например, для оптимизации многоступенчатых ракет [1].
Многоступенчатая ракета — это аппарат, в котором части конструкции отделяются во время полета, придавая оставшейся части ракеты дополнительную скорость. Трёхступенчатая ракета схематически показана на рисунке.
![](https://habrastorage.org/webt/lw/fp/kc/lwfpkcmed2d0lhdigvold_76qgy.png)
По мере движения ракеты, ступени отделяются до тех пор, пока не останется главная часть ракеты, несущая полезную нагрузку. Задача оптимизации ракеты состоит в таком распределении веса по ступеням, при котором определенная целевая функция достигает максимального либо минимального значения.
Мы рассмотрим две задачи в предположении, что коэффициент
![](https://habrastorage.org/webt/ka/eb/uv/kaebuvntumpsdwf4-ff8ivuiubu.png)
Характеристики многоступенчатой ракеты можно описать двумя уравнениями. Первое уравнение для коэффициента полезной нагрузки ракеты:
![](https://habrastorage.org/webt/r2/lv/ch/r2lvchdx5w7_ngo-mhjufelbg7e.png)
где: W1– полезный вес ракеты ;WN –начальный вес ракеты до отделения ступеней.
Второе уравнение для характеристики полезной нагрузки многоступенчатой ракеты, которое можно записать как выражение для е? идеальной конечной скорости (формула К.Э.Циолковского):
![](https://habrastorage.org/webt/ou/e9/by/oue9byv7kef7v0swkxve9objmmk.png)
где:
![](https://habrastorage.org/webt/zr/ua/ux/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
Представим формулу для вычисления заданной скорости в виде:
![](https://habrastorage.org/webt/6x/u4/c2/6xu4c2isnnhz3ztbcf6tsih8-rk.png)
где
![](https://habrastorage.org/webt/gb/xa/bc/gbxabc9nu2gnf18vfdlrfdobcdq.png)
Отсюда следует, что:
![](https://habrastorage.org/webt/7u/tm/zg/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()
Получим:
![](https://habrastorage.org/webt/ed/zi/0w/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
![](https://habrastorage.org/webt/ly/do/km/lydokm4hs7wpnocgy0dvvddff78.png)
Вес ракеты распределен так (4000.0; 1967.0; 967.9; 476.1), что соответствует нарастанию скорости равными долями (после отделения каждой ступени) в среднем на 460 м/с.
Используем ту же формулу для конечной скорости ракеты, что и в предыдущей задаче:
![](https://habrastorage.org/webt/j7/7p/ej/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()
Получим контурный график функции цели с отрицательным знаком для смены минимума на максимум:
![](https://habrastorage.org/webt/br/bh/2p/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
![](https://habrastorage.org/webt/tl/kl/2k/tlkl2kxzpai1tibd24awmoywu0c.png)
Веса ракеты распределяются так :4000;1473.625;542.889;200, что соответствует нарастанию скорости равными долями (после отделения каждой ступени) на 526.9 м/с.
Выводы
В публикации рассмотрены особенности применения библиотеки scipy. optimize при решении задач оптимизации на примере оптимизации многоступенчатой ракеты. Полученные результаты в части неизменности начальной точки поиска экстремума для переменных и формирования функции цели в привычном пользователю виде будут способствовать расширению области применения свободно распространяемого языка программирования Python.
Всем спасибо за внимание! И главное, чтобы многоступенчатая ракета с не полезной для нас нагрузкой НИКОГДА к нам не прилетела.
Ссылки
Комментарии (8)
nomadmoon
10.12.2017 18:41Всегда думал что части от ракеты отваливаются не для того чтобы увеличить скорость оставшейся части а чтобы уменьшить её массу. Соответственно ступенчатые ракеты имеет смысл строить для взлёта из гравитационного колодца, а сферическая ракета в вакууме она, хмм, сферическая.
encyclopedist
10.12.2017 20:45Даже без гравитации нет никакого смысла тратить энергию на разгон бесполезных уже баков.
И да, преодоление гравитации — очень маленькая часть выведения даже на НОО.
Eklykti
10.12.2017 21:44Сферическая ракета в вакууме будет работать только в том случае, если у неё ЭРД с высоким УИ, но в таком случае разгоняться она будет ну оооооочень долго, по несколько месяцев.
arheops
11.12.2017 13:28Сферическая ракета в вакууме по вашей схеме будет работать только если у нее нелинейный двигатель. А именно, если тяга двигателя прямо пропорциональна массе в данной точке(такой себе хитрый антигравитатор).
В остальных случаях вам понадобится больше энергии на разворот и/или разгон большей массы.
Daniil1979
10.12.2017 21:21За пример решения математической задачи на Питоне спасибо. И за пример использования SciPy отдельное спасибо.
KiloLeo
10.12.2017 22:21В такой постановке задача вполне решается аналитически, не видно смысла решать её численно. Построили целевую функцию — взяли производную — нашли экстремумы. Численное решение оправданно если задача аналитически не решается.
datacompboy
Родовая болячка таких задач — постановка задач без объяснения «что это значит».
«минимизировать коэффициент полезной нагрузки» — это что означает «на пальцах»?
«максимизировать идеальную конечную скорость» — это как? идеальная конечная скорость она ж идеальна, как она может быть больше или меньше?