Untitled
unknown
plain_text
2 years ago
24 kB
3
Indexable
#/////////////////////////////////////Используемые библиотеки///////////////////////////////////////////////////////////
import numpy as np
import matplotlib.pyplot as plt
#///////////////////////////////////////Используемые функции////////////////////////////////////////////////////////////
FINAL = 0
# Уравнения движения средства выведения (СВ) в стартовой системе координат (Oxy):
def dVx(t, x, y, P, m, ϑ1_t, ϑ2, t_1, t_2):
return P * np.cos(ϑ(t, ϑ1_t, ϑ2, t_1, t_2))/m + gx(x, y)
def dVy(t, x, y, P, m, ϑ1_t, ϑ2, t_1, t_2):
return P * np.sin(ϑ(t, ϑ1_t, ϑ2, t_1, t_2))/m + gy(x, y)
def dx(Vx):
return Vx
def dy(Vy):
return Vy
def dm(P):
return - (P / W_UCT)
# Проекции ускорения от гравитационного поля Луны на соответствующие оси стартовой СК
def gx(x, y):
return - (µ_m * x)/((x**2 + (R_m + y)**2)**(3/2))
def gy(x, y):
return - (µ_m * (R_m + y))/((x**2 + (R_m + y)**2)**(3/2))
# Метод Рунге-Кутта 4-го порядка
def RK4(t, x, y, Vx, Vy, P, m, ϑ1_t, ϑ2, dt, t_1, t_2):
# K1
'''if (t_1 < t <= t_2):
P = 0'''
Vx_1 = dVx(t, x, y, P, m, ϑ1_t, ϑ2, t_1, t_2) * dt
Vy_1 = dVy(t, x, y, P, m, ϑ1_t, ϑ2, t_1, t_2) * dt
x_1 = dx(Vx) * dt
y_1 = dy(Vy) * dt
m_1 = dm(P) * dt
# K2
'''if (t == t_1):
P = 0
if (t == t_2):
P = P_dano'''
Vx_2 = dVx(t + dt / 2, x + x_1 / 2, y + y_1 / 2, P, m + m_1 / 2, ϑ1_t, ϑ2, t_1, t_2) * dt
Vy_2 = dVy(t + dt / 2, x + x_1 / 2, y + y_1 / 2, P, m + m_1 / 2, ϑ1_t, ϑ2, t_1, t_2) * dt
x_2 = dx(Vx + Vx_1 / 2) * dt
y_2 = dy(Vy + Vy_1 / 2) * dt
m_2 = dm(P) * dt
# K3
Vx_3 = dVx(t + dt / 2, x + x_2 / 2, y + y_2 / 2, P, m + m_2 / 2, ϑ1_t, ϑ2, t_1, t_2) * dt
Vy_3 = dVy(t + dt / 2, x + x_2 / 2, y + y_2 / 2, P, m + m_2 / 2, ϑ1_t, ϑ2, t_1, t_2) * dt
x_3 = dx(Vx + Vx_2 / 2) * dt
y_3 = dy(Vy + Vy_2 / 2) * dt
m_3 = dm(P) * dt
# K4
Vx_4 = dVx(t + dt, x + x_3, y + y_3, P, m + m_3, ϑ1_t, ϑ2, t_1, t_2) * dt
Vy_4 = dVy(t + dt, x + x_3, y + y_3, P, m + m_3, ϑ1_t, ϑ2, t_1, t_2) * dt
x_4 = dx(Vx + Vx_3) * dt
y_4 = dy(Vy + Vy_3) * dt
m_4 = dm(P) * dt
# Нач. значение + усредненное значение приращений у
Vx += (Vx_1 + Vx_2 * 2 + Vx_3 * 2 + Vx_4) / 6
Vy += (Vy_1 + Vy_2 * 2 + Vy_3 * 2 + Vy_4) / 6
x += (x_1 + x_2 * 2 + x_3 * 2 + x_4) / 6
y += (y_1 + y_2 * 2 + y_3 * 2 + y_4) / 6
m += (m_1 + m_2 * 2 + m_3 * 2 + m_4) / 6
return Vx, Vy, x, y, m
# Функция интегрирования
def intDY(P, ϑ1_t, ϑ2, t_1, t_2, h, dt):
t, x, y, Vx, Vy = [0, 0, 0, 0, 0]
m = m0
x_List = []
y_List = []
# Первый активный участок
# Вертикальный
while t < t_B:
Vx, Vy, x, y, m = RK4(t, x, y, Vx, Vy, P, m, ϑ1_t, ϑ2, dt, t_1, t_2)
t += dt
if FINAL == 1:
x_List.append(x)
y_List.append(y)
if round(t,4)==t_B:
if FINAL == 1:
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', t)
# График
if FINAL == 1:
grafik(x_List, y_List, 'grey')
# Не вертикальный
while t < t_1:
if round(t, 1) == round(t_1 - (t_1 % 0.1), 1):
Vx, Vy, x, y, m = RK4(t, x, y, Vx, Vy, P, m, ϑ1_t, ϑ2, t_1-t, t_1, t_2)
if FINAL == 1:
x_List.append(x)
y_List.append(y)
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', round(t_1, 3))
# Параметры на точке 360.456
# Находим параметры для 360.5
Vx, Vy, x, y, m = RK4(t, x, y, Vx, Vy, 0, m, ϑ1_t, ϑ2, Dt-(t_1-t), t_1, t_2)
if FINAL == 1:
x_List.append(x)
y_List.append(y)
t += Dt
if FINAL == 1:
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', t)
else:
Vx, Vy, x, y, m = RK4(t, x, y, Vx, Vy, P, m, ϑ1_t, ϑ2, dt, t_1, t_2)
if FINAL == 1:
x_List.append(x)
y_List.append(y)
t += dt
if round(t, 1) == round(t_1 - (t_1 % 0.1), 1):
if FINAL == 1:
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', t)
# График
if FINAL == 1:
grafik(x_List, y_List, 'orangered')
# Второй пасивный участок
while t < t_2:
if round(t, 1) == round(t_2 - (t_2 % 0.1), 1):
Vx, Vy, x, y, m = RK4(t, x, y, Vx, Vy, 0, m, ϑ1_t, ϑ2, t_2-t, t_1, t_2)
if FINAL == 1:
x_List.append(x)
y_List.append(y)
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', round(t_2, 3))
# Находим параметры для точки нового участка
Vx, Vy, x, y, m = RK4(t, x, y, Vx, Vy, P, m, ϑ1_t, ϑ2, Dt-(t_2-t), t_1, t_2)
if FINAL == 1:
x_List.append(x)
y_List.append(y)
t += Dt
if FINAL == 1:
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', t)
else:
Vx, Vy, x, y, m = RK4(t, x, y, Vx, Vy, 0, m, ϑ1_t, ϑ2, dt, t_1, t_2)
if FINAL == 1:
x_List.append(x)
y_List.append(y)
t += dt
if round(t, 1) == round(t_2 - (t_2 % 0.1), 1):
if FINAL == 1:
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', t)
# График
if FINAL == 1:
grafik(x_List, y_List, 'lightblue')
V_UCL = np.sqrt(µ_m / (R_m + h)) # необходимая скорость
if np.sqrt(Vx ** 2 + Vy ** 2)>V_UCL:
print('ошибка')
while abs(np.sqrt(Vx ** 2 + Vy ** 2) - (V_UCL)) > (10 ** (-8)):
if abs(np.sqrt(Vx ** 2 + Vy ** 2) - (V_UCL)) <= (10 ** (-8)):
break
Vx_prior = Vx
Vy_prior = Vy
x_prior = x
y_prior = y
m_prior = m
Vx, Vy, x, y, m = RK4(t, x, y, Vx, Vy, P, m, ϑ1_t, ϑ2, dt, t_1, t_2)
if np.sqrt(Vx ** 2 + Vy ** 2) > (V_UCL):
Vx = Vx_prior
Vy = Vy_prior
x = x_prior
y = y_prior
m = m_prior
dt /= 10
# print(np.sqrt(Vx**2 + Vy**2))
else:
if FINAL == 1:
x_List.append(x)
y_List.append(y)
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', t)
t += dt
'''if round(m, ) < 1465:
Vx = Vx_prior
Vy = Vy_prior
x = x_prior
y = y_prior
m = m_prior
t -= dt
break'''
if FINAL == 1:
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', t)
print(np.sqrt(Vx ** 2 + Vy ** 2))
print(V_UCL)
# График
if FINAL == 1:
grafik(x_List, y_List, 'palegreen')
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(0, 900000)
plt.ylim(0, 350000)
plt.title('Trajectory')
# plt.show()
return Vx, Vy, x, y, m, t
# Графики частей траекторий
def grafik(x_List, y_List, color ='red'):
plt.plot(x_List[:-1], y_List[:-1], color=color)
plt.plot(x_List[-1:], y_List[-1:], marker='o', markerfacecolor='blue', markersize=8)
x_List.clear()
y_List.clear()
# Интегрирование
# Функция краевой задачи
def r_and_ϑ(Vx, Vy, x, y):
V = np.sqrt(Vx ** 2 + Vy ** 2)
r = np.sqrt((R_m + y) ** 2 + x ** 2)
yϑ = np.arcsin((Vx * x + Vy * (R_m + y)) / (r * V))
return r, yϑ
def CHASTNAYA_PROIZVOD(f2, f1, delta):
return (f2 - f1) / (2 * delta)
def KRAEVAYA(P, ϑ1_t, ϑ2, t_1, t_2, h, dt):
# шаг
d_ϑ1_t = 10**(-8)
d_ϑ2 = 10**(-5)
r_otl = 1
ϑ_otl = 1
# Итерационная переменная
U_i_new = np.array([[ϑ1_t], [ϑ2]])
while ((r_otl / (10**(-3)))**2 + (ϑ_otl / (10**(-5)))**2)**0.5 > 1:
if FINAL == 1:
tz = 0
#print("U_i_new", U_i_new)
# Изменение углов
U_i = U_i_new
ϑ1_t = U_i[0][0]
ϑ2 = U_i[1][0]
# Для текущих значений
Vx, Vy, x, y, m, t = intDY(P, ϑ1_t, ϑ2, t_1, t_2, h, dt)
r_tek, ϑ_tek = r_and_ϑ(Vx, Vy, x, y)
# Для производной
# r и ϑ при ϑ1_t - d_ϑ1_t
Vx, Vy, x, y, m, t = intDY(P, ϑ1_t - d_ϑ1_t, ϑ2, t_1, t_2, h, dt)
r_ϑ1_t_min_d_ϑ1_t, ϑ_ϑ1_t_min_d_ϑ1_t = r_and_ϑ(Vx, Vy, x, y)
# r и ϑ при ϑ1_t + d_ϑ1_t
Vx, Vy, x, y, m, t = intDY(P, ϑ1_t + d_ϑ1_t, ϑ2, t_1, t_2, h, dt)
r_ϑ1_t_pls_d_ϑ1_t, ϑ_ϑ1_t_pls_d_ϑ1_t = r_and_ϑ(Vx, Vy, x, y)
# r и ϑ при ϑ2 - d_ϑ2
Vx, Vy, x, y, m, t = intDY(P, ϑ1_t, ϑ2 - d_ϑ2, t_1, t_2, h, dt)
r_ϑ2_min_d_ϑ2, ϑ_ϑ2_min_d_ϑ2 = r_and_ϑ(Vx, Vy, x, y)
# r и ϑ при ϑ2 + d_ϑ2
Vx, Vy, x, y, m, t = intDY(P, ϑ1_t, ϑ2 + d_ϑ2, t_1, t_2, h, dt)
r_ϑ2_pls_d_ϑ2, ϑ_ϑ2_pls_d_ϑ2 = r_and_ϑ(Vx, Vy, x, y)
# Производные
# dr при изм. ϑ1_t
dr_ϑ1_t = CHASTNAYA_PROIZVOD(r_ϑ1_t_pls_d_ϑ1_t, r_ϑ1_t_min_d_ϑ1_t, d_ϑ1_t)
# dr при изм. ϑ2
dr_ϑ2 = CHASTNAYA_PROIZVOD(r_ϑ2_pls_d_ϑ2, r_ϑ2_min_d_ϑ2, d_ϑ2)
# dϑ при изм. ϑ1_t
dϑ_ϑ1_t = CHASTNAYA_PROIZVOD(ϑ_ϑ1_t_pls_d_ϑ1_t, ϑ_ϑ1_t_min_d_ϑ1_t, d_ϑ1_t)
# dϑ при изм. ϑ2
dϑ_ϑ2 = CHASTNAYA_PROIZVOD(ϑ_ϑ2_pls_d_ϑ2, ϑ_ϑ2_min_d_ϑ2, d_ϑ2)
# Считаем Якобиан
J = np.array([[dr_ϑ1_t, dr_ϑ2], [dϑ_ϑ1_t, dϑ_ϑ2]])
# Находим обратную
J_1 = np.linalg.inv(J)
# Матрица невязок
r_otl = r_tek - (R_m + h)
ϑ_otl = ϑ_tek
F_U = - np.array([[r_otl],[ϑ_otl]])
if FINAL == 1:
tz = 0
#print("F_U", F_U)
# Матрица искомых ϑ1_t и ϑ2
U_i_new = U_i + J_1 @ F_U
return U_i_new
def iteration_formula_Levevi(X_i, H_i, alpha_i, I, grad_g):
return X_i - np.linalg.inv(H_i + alpha_i * I) @ grad_g
#//////////////////////////////////////////Исходные данные//////////////////////////////////////////////////////////////
# N = 11
Pj = [0,8250,9920]
W_UCT = 3520
m0 = 3015
Dm_Tmax = 1550
m_KOHCT = 615
t_B = 14
Dt = 0.1
h_UCLj = [0,125000,170000,170000,273000]
#/////////////////////////////////////////////Допущения/////////////////////////////////////////////////////////////////
R_m = 1738000 # радиус Луны
µ_m = 4.903 * (10**12) # гравитационная параметр Луны
#/////////// Параметры закона управления ориентацией вектора тяги на АУТ и параметры циклограммы работы ДУ//////////////
# Производная и угол тангажа для управления вектором тяги
ϑ1_t = -0.003 # radians
ϑ2 = -0.05 # radians
# Время для участков
t_1 = 350.364
t_2 = 600.958
# Закон изменения угла тангажа на активных участках траектории
def ϑ(t, ϑ1_t, ϑ2, t_1, t_2):
if t <= t_B:
ϑ_Y4ACTKA = np.pi / 2
elif (t_B < t <= t_1):
ϑ_Y4ACTKA = np.pi/2 + ϑ1_t * (t - t_B)
elif (t > t_2):
ϑ_Y4ACTKA = ϑ2
else:
ϑ_Y4ACTKA = 0
return ϑ_Y4ACTKA
def OPTIMIZACIYA(accuracy = 0.01, t_1 = t_1, t_2 = t_2):
# Номер итерации
iterations = 0
print('Итерация №', iterations)
C_1 = 0.5
C_2 = 1/C_1
# Матрица для рекурентной формулы
I = np.array([1, 0, 0, 1])
I.shape = (2, 2)
# Величина шагов
dt1 = 0.01
dt2 = 0.01
epsilon = accuracy
alpha_i = 0.001
print('Считаем краевую без оптимизации')
# Краевая без оптимизации
U_i_new = KRAEVAYA(P, ϑ1_t, ϑ2, t_1, t_2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_new[0][0], U_i_new[1][0], t_1, t_2, h, Dt)
# Функция g
g_new = m0 - m
print('30%')
print('Считаем градиент')
# Считаем градиент
# t1 - dt1
U_i_grad_11 = KRAEVAYA(P, ϑ1_t, ϑ2, t_1 - dt1, t_2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_grad_11[0][0], U_i_grad_11[1][0], t_1 - dt1, t_2, h, Dt)
g_grad_11 = m0 - m
# t1 + dt1
U_i_grad_12 = KRAEVAYA(P, ϑ1_t, ϑ2, t_1 + dt1, t_2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_grad_12[0][0], U_i_grad_12[1][0], t_1 + dt1, t_2, h, Dt)
g_grad_12 = m0 - m
# Частная производная
g_grad_t1 = CHASTNAYA_PROIZVOD(g_grad_12, g_grad_11, dt1)
print('Посчитали частную производную для t_1 50%')
# t2 - dt2
U_i_grad_21 = KRAEVAYA(P, ϑ1_t, ϑ2, t_1, t_2 - dt2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_grad_21[0][0], U_i_grad_21[1][0], t_1, t_2 - dt2, h, Dt)
g_grad_21 = m0 - m
# t2 + dt2
U_i_grad_22 = KRAEVAYA(P, ϑ1_t, ϑ2, t_1, t_2 + dt2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_grad_22[0][0], U_i_grad_22[1][0], t_1, t_2 + dt2, h, Dt)
g_grad_22 = m0 - m
# Частная производная
g_grad_t2 = CHASTNAYA_PROIZVOD(g_grad_22, g_grad_21, dt2)
print('Посчитали частную производную для t_2 99%')
# Градиент
grad_g = np.array([g_grad_t1, g_grad_t2])
grad_g.shape = (2, 1)
print('Нач. Градиент g = ', grad_g)
print('Посчитали нач. градиент 100%')
# Чтобы войти в цикл
X_new = np.array([t_1, t_2])
X_new.shape = (2, 1)
# Пока модуль градиента больше точности
while abs(np.sqrt(grad_g[0] ** 2 + grad_g[1] ** 2)) >= epsilon:
iterations += 1
print('Итерация №', iterations)
# Переобозначения t_1 и t_2
t_1_p = X_new[0][0]
t_2_p = X_new[1][0]
g = g_new
# Смешанная производная
print('Считаем смешанные производные')
# значение в 1 точке
U_i_1 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], t_1_p - dt1, t_2_p - 2 * dt2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_1[0][0], U_i_1[1][0], t_1_p - dt1, t_2_p - 2 * dt2, h, Dt)
g_1 = m0 - m
# значение во 2 точке
U_i_2 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], t_1_p - dt1, t_2_p + 2 * dt2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_2[0][0], U_i_2[1][0], t_1_p - dt1, t_2_p + 2 * dt2, h, Dt)
g_2 = m0 - m
# значение в 3 точке
U_i_3 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], t_1_p + dt1, t_2_p - 2 * dt2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_3[0][0], U_i_3[1][0], t_1_p + dt1, t_2_p - 2 * dt2, h, Dt)
g_3 = m0 - m
# значение в 4 точке
U_i_4 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], t_1_p + dt1, t_2_p + 2 * dt2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_4[0][0], U_i_4[1][0], t_1_p + dt1, t_2_p + 2 * dt2, h, Dt)
g_4 = m0 - m
# Находим 1-ые производные
f1_dt1_der = CHASTNAYA_PROIZVOD(g_3, g_1, dt2)
f2_dt1_der = CHASTNAYA_PROIZVOD(g_4, g_2, dt2)
print('Посчитали 1-вые производные 10%')
# Находим 2 производную
f_t1_t2_der = CHASTNAYA_PROIZVOD(f2_dt1_der, f1_dt1_der, 2 * dt1)
# 2 3 элемент Гессиана
print('Посчитали 2-ю производную 20%')
# Полная производная
print('Считаем полную производную')
# Для оси t_1
# Значение в 1 точке
U_i_11 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], t_1_p - dt1 - dt2, t_2_p, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_11[0][0], U_i_11[1][0], t_1_p - dt1 - dt2, t_2_p, h, Dt)
g_11 = m0 - m
# Значение во 2 точке
U_i_12 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], t_1_p - dt1 + dt2, t_2_p, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_12[0][0], U_i_12[1][0], t_1_p - dt1 + dt2, t_2_p, h, Dt)
g_12 = m0 - m
# Значение в 3 точке
U_i_13 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], t_1_p + dt1 - dt2, t_2_p, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_13[0][0], U_i_13[1][0], t_1_p + dt1 - dt2, t_2_p, h, Dt)
g_13 = m0 - m
# Значение в 4 точке
U_i_14 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], t_1_p + dt1 + dt2, t_2_p, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_14[0][0], U_i_14[1][0], t_1_p + dt1 + dt2, t_2_p, h, Dt)
g_14 = m0 - m
# Находим 1-ые производные
f11_dt1_der = CHASTNAYA_PROIZVOD(g_11, g_12, dt2)
f12_dt1_der = CHASTNAYA_PROIZVOD(g_13, g_14, dt2)
print('Посчитали 1-вые производные для t_1 30%')
# Находим 2-ую производную
f_t1_t1_der = CHASTNAYA_PROIZVOD(f11_dt1_der, f12_dt1_der, dt1)
print('Посчитали 2-ю производную для t_1 35%')
# Для оси t_2
# Значение в 1 точке
U_i_21 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], t_1_p, t_2_p - dt2 - dt1, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_21[0][0], U_i_21[1][0], t_1_p, t_2_p - dt2 - dt1, h, Dt)
g_21 = m0 - m
# Значение во 2 точке
U_i_22 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], t_1_p, t_2_p - dt1 + dt2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_22[0][0], U_i_22[1][0], t_1_p, t_2_p - dt1 + dt2, h, Dt)
g_22 = m0 - m
# Значение в 3 точке
U_i_23 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], t_1_p, t_2_p + dt1 - dt2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_23[0][0], U_i_23[1][0], t_1_p, t_2_p + dt1 - dt2, h, Dt)
g_23 = m0 - m
# Значение в 4 точке
U_i_24 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], t_1_p, t_2_p + dt2 + dt1, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_24[0][0], U_i_24[1][0], t_1_p, t_2_p + dt2 + dt1, h, Dt)
g_24 = m0 - m
# Находим 1-ые производные
f21_dt2_der = CHASTNAYA_PROIZVOD(g_21, g_22, dt1)
f22_dt2_der = CHASTNAYA_PROIZVOD(g_23, g_24, dt1)
print('Посчитали 1-вые производные для t_2 40%')
# Находим 2-ую производную
f_t2_t2_der = CHASTNAYA_PROIZVOD(f21_dt2_der, f22_dt2_der, dt2)
print('Посчитали 2-ю производную для t_2 45%')
# Чем меньше альфа, тем быстрее шагаем
# Гессиан
#H_i = np.array([f_t1_t1_der, f_t1_t2_der, f_t2_t1_der, f_t2_t2_der])
H_i = np.array([f_t1_t1_der, f_t1_t2_der, f_t1_t2_der, f_t2_t2_der])
H_i.shape = (2, 2)
print('Гессиан =', H_i)
print('Посчитали Гессиан 50%')
C_2_flag = False
# Пока новая g >= g
while g_new >= g:
X = X_new
print('Шаг к g_min')
# Новая матрица времён
X_new = iteration_formula_Levevi(X_new, H_i, alpha_i, I, grad_g)
U_i_new = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], X_new[0][0], X_new[1][0], h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_new[0][0], U_i_new[1][0], X_new[0][0], X_new[1][0], h, Dt)
g_new = m0 - m
# Проверяем условие
if g_new < g:
if C_2_flag == False:
alpha_i *= C_1
print('1', alpha_i)
break
# Если перешагнули
alpha_i *= C_2
X_new = X
print('2', alpha_i)
C_2_flag = True
print('X_new = ', X_new)
print('U_i_new = ', U_i_new)
print('g_new = ', g_new)
print('Считаем новый градиент')
# Считаем новый градиент
# t1 - dt1
U_i_recalc_grad_11 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], X_new[0][0] - dt1, X_new[1][0], h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_recalc_grad_11[0][0], U_i_recalc_grad_11[1][0], X_new[0][0] - dt1, X_new[1][0], h, Dt)
g_recalc_grad_11 = m0 - m
# t1 + dt1
U_i_recalc_grad_12 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], X_new[0][0] + dt1, X_new[1][0], h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_recalc_grad_12[0][0], U_i_recalc_grad_12[1][0], X_new[0][0] + dt1, X_new[1][0], h, Dt)
g_recalc_grad_12 = m0 - m
# Частная производная
g_recalc_grad_t1 = CHASTNAYA_PROIZVOD(g_recalc_grad_12, g_recalc_grad_11, dt1)
print('Посчитали частную производную для t_1 80%')
# t2 - dt2
U_i_recalc_grad_21 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], X_new[0][0], X_new[1][0] - dt2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_recalc_grad_21[0][0], U_i_recalc_grad_21[1][0], X_new[0][0], X_new[1][0] - dt2, h, Dt)
g_recalc_grad_21 = m0 - m
# t2 + dt2
U_i_recalc_grad_22 = KRAEVAYA(P, U_i_new[0][0], U_i_new[1][0], X_new[0][0], X_new[1][0] + dt2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i_recalc_grad_22[0][0], U_i_recalc_grad_22[1][0], X_new[0][0], X_new[1][0] + dt2, h, Dt)
g_recalc_grad_22 = m0 - m
# Частная производная
g_recalc_grad_t2 = CHASTNAYA_PROIZVOD(g_recalc_grad_22, g_recalc_grad_21, dt2)
print('Посчитали частную производную для t_2 90%')
grad_g = np.array([g_recalc_grad_t1, g_recalc_grad_t2])
grad_g.shape = (2, 1)
print('Посчитали нов. градиент 95%', grad_g)
print('mod(grad_g) =', abs(np.sqrt(grad_g[0] ** 2 + grad_g[1] ** 2)))
print('Итерация №',str(iterations),'завершена 100%')
Vx, Vy, x, y, m, t = intDY(P, U_i_new[0][0], U_i_new[1][0], X_new[0][0],X_new[1][0], h, Dt)
return X_new, U_i_new
#/////////////////////////////Интегрирование системы дифференциальных уравнений/////////////////////////////////////////
P = Pj[2]
P_dano = P
h = h_UCLj[1]
# intDY(P, ϑ1_t, ϑ2, t_1, t_2, h, Dt)
#//////////////////////////////////////////////Краевая задача///////////////////////////////////////////////////////////
'''FINAL = 1
U_i = KRAEVAYA(P, ϑ1_t, ϑ2, t_1, t_2, h, Dt)
Vx, Vy, x, y, m, t = intDY(P, U_i[0][0], U_i[1][0], t_1, t_2, h, Dt)'''
#////////////////////////////////////////////Задача оптимизации/////////////////////////////////////////////////////////
X_i, U_i = OPTIMIZACIYA()
FINAL = 1
Vx, Vy, x, y, m, t = intDY(P, U_i[0][0], U_i[1][0], X_i[0][0], X_i[1][0], h, Dt)
plt.show()Editor is loading...
Leave a Comment