laba5 mo

# численные методы поиска условного эсктремума
# стр. 105
import numpy as np
import math as m

import scipy.optimize
from sympy import *
from scipy.optimize import linprog

import Lab4

x1 = Symbol('x1')
x2 = Symbol('x2')
y1 = Symbol('y1')
y2 = Symbol('y2')
z = Symbol('z')
betha = Symbol('betha')


def solve_system(g_functions, x_point):
    result = []
    for i in range(0, len(g_functions)):
        g = lambdify((x1, x2), g_functions[i], 'numpy')
        if g(x_point[0], x_point[1]) <= 0:
            result.append(i)
    return result


def get_alpha_value(g_functions: list[Function], f: Function, x_point: list[float], y_vector: list[float]):
    f_func = lambdify((x1, x2), f, 'numpy')
    f_betha = f_func(x_[0] + y[0] * betha, x_[1] + y[1] * betha)
    f_betha = lambdify(betha, f_betha, 'numpy')
    betha0 = scipy.optimize.fmin(f_betha, 0)
    bethas = [betha0]
    # Значений для betha и g функций с ней
    betha_functions = []
    betha_variables = []
    i = 1
    for g in g_functions:
        # Переменная бета, т.е. бета1 бета2 бета3 и т.д.
        betha_variable = Symbol('betha' + string(i))
        g_func = lambdify((x1, x2), g, 'numpy')
        g_func = g_func(x_point[0] + y_vector[0] * betha_variable, x_point[1] + y_vector[1] * betha_variable)
        betha_functions.append(g_func)
        betha_variables.append(betha_variable)
        i = i + 1
    result = nonlinsolve(betha_functions, betha_variables)
    return min(result)


def possible_directions_method(f: Function, g_functions, x0, epsilon):
    # Шаг 1
    x_ = np.array(x0)
    k = 0
    while k < 1000:
        # Шаг 2
        f_gradient = Lab4.gradient(f, x1, x2, x_)
        if np.linalg.norm(f_gradient) <= epsilon:

            # получение решения #
            return True, x_, None

        else:
            # Шаг 3
            i_set = solve_system(g_functions, x_)
            # Шаг 4
            if len(i_set) == 0:
                yk = -1 * Lab4.gradient(f, x1, x2, x_)
            else:
                a = []
                b = [0]
                c = [0, 0, 1]  # значения коэффициентов функции для y1 y2 z
                row = Lab4.gradient(f, x1, x2, x_)
                row.append(-1)  # задание значений для y1, y2 и -z
                a.append(row)
                for i in i_set:
                    if i != 0:
                        row = Lab4.gradient(g_functions[i], x1, x2, x_)  # значения для y1 y2
                        row.append(-1)  # значения для -z
                        a.append(row)  # формирования коэффициентов для неравенств с j-ми функциями
                        b.append(0)  # формирование ограничений для j-ых функций
                y1_bounds = [-1, 1]  # ограничения для y1
                y2_bounds = [-1, 1]  # ограничения для y2
                z_bounds = [None, None]
                bounds = [y1_bounds, y2_bounds, z_bounds]
                linprog_result = linprog(c, A_ub=a, b_ub=b, bounds=bounds, method='simplex', options={"disp": True})
                print('Результаты решения симплекс метода:\n', linprog_result)
                # шаг 5
                yk = [linprog_result.x[0], linprog_result.x[1]]
                if yk == [0, 0] or linprog_result.x[2] <= epsilon:
                    # получение решения #
                    return True, x_, yk
                else:
                    yk = [linprog_result.x[0], linprog_result.x[1]]
                    # получение минимума из всех бета, что = альфа #
                    alpha_k = get_alpha_value(g_functions, f, x_, yk)

                    # шаг 6
                    # итерация #
                    x_ = x_ + alpha_k * np.array(yk)
                    k = k + 1


def main():
    x0 = [0, 0.95]
    epsilon = 0.03
    f = (x1 - 4) ** 2 + (x2 - 5) ** 2
    g = [x1 + x2 - 1, x1, x2]
    solved, x_result, y_result = possible_directions_method(f, g, x0, epsilon)
    if solved:
        print("Задача была решена, точка x = ", x_result)
        if y_result is not None:
            print("Вектор y = ", y_result)
    else:
        print("Задача не была решена")


if __name__ == "__main__":
    main()