Untitled

mail@pastecode.io avatar
unknown
plain_text
7 months ago
24 kB
0
Indexable
Never
#/////////////////////////////////////Используемые библиотеки///////////////////////////////////////////////////////////
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):
    yϑ = ϑ(t, ϑ1_t, ϑ2, t_1, t_2)
    dVx = P * np.cos(yϑ)/m + gx(x, y)
    return dVx
def dVy(t, x, y, P, m, ϑ1_t, ϑ2, t_1, t_2):
    yϑ = ϑ(t, ϑ1_t, ϑ2, t_1, t_2)
    dVy = P * np.sin(yϑ)/m + gy(x, y)
    return dVy
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):

    delta_t_1 = round(t_1 % 0.1, 3)
    t_1_to4ny = round(t_1 - delta_t_1, 2)
    delta_t_2 = round(t_2 % 0.1, 3)
    t_2_to4ny = round(t_2 - delta_t_2, 2)

    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:
                tz=0
                #print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', t)
    # График
    if FINAL == 1:
        grafik(x_List, y_List, 'grey')
    # Не вертикальный
    while t < t_1:
        if round(t, 2) == t_1_to4ny:
            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, P, 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:
                tz = 0
                #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, 2) == t_1_to4ny:
                if FINAL == 1:
                    tz = 0
                    #print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', t)
    # График
    if FINAL == 1:
        grafik(x_List, y_List, 'orangered')
    # Второй пасивный участок
    while t < t_2:
        if round(t, 2) == t_2_to4ny:
            Vx, Vy, x, y, m = RK4(t, x, y, Vx, Vy, P, 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_1-t), t_1, t_2)
            if FINAL == 1:
                x_List.append(x)
                y_List.append(y)
            t += Dt
            if FINAL == 1:
                tz = 0
                #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, 2) == t_2_to4ny:
                if FINAL == 1:
                    tz = 0
                    #print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', t)
    # График
    if FINAL == 1:
        grafik(x_List, y_List, 'lightblue')
    V_UCL = np.sqrt(µ_m / (R_m + h))  # необходимая скорость
    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:
        tz = 0
        #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 = ϑ1_t / 10
    d_ϑ2 = ϑ2 / 10

    r_otl = 1
    ϑ_otl = 1

    # Итерационная переменная
    U_i_new = np.array([[ϑ1_t], [ϑ2]])

    while abs(round(ϑ_otl, 3)) > 10**(-3) or abs(round(r_otl, 3)) > 10**(-3):
        if FINAL == 1:
            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:
            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.04 # 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_nach = 350.364
    t_2_nach = 600.958
    # Номер итерации
    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_nach, t_2_nach, h, Dt)
    Vx, Vy, x, y, m, t = intDY(P, U_i_new[0][0], U_i_new[1][0], t_1_nach, t_2_nach, h, Dt)
    # Функция g
    g_new = m0 - m
    print('30%')
    print('Считаем градиент')
    # Считаем градиент
    # t1 - dt1
    U_i_grad_11 = KRAEVAYA(P, ϑ1_t, ϑ2, t_1_nach - dt1, t_2_nach, h, Dt)
    Vx, Vy, x, y, m, t = intDY(P, U_i_grad_11[0][0], U_i_grad_11[1][0], t_1_nach - dt1, t_2_nach, h, Dt)
    g_grad_11 = m0 - m
    # t1 + dt1
    U_i_grad_12 = KRAEVAYA(P, ϑ1_t, ϑ2, t_1_nach + dt1, t_2_nach, h, Dt)
    Vx, Vy, x, y, m, t = intDY(P, U_i_grad_12[0][0], U_i_grad_12[1][0], t_1_nach + dt1, t_2_nach, 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_nach, t_2_nach - 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_nach, t_2_nach - dt2, h, Dt)
    g_grad_21 = m0 - m
    # t2 + dt2
    U_i_grad_22 = KRAEVAYA(P, ϑ1_t, ϑ2, t_1_nach, t_2_nach + 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_nach, t_2_nach + 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_nach, t_2_nach])
    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, 2 * dt2)
        f2_dt1_der = CHASTNAYA_PROIZVOD(g_4, g_2, 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[1]
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()
Leave a Comment