Lab5

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

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


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


def fd(func, x_to_diff, x_to_lambdify):
    fdx = func.diff(x_to_diff)
    return lambdify(x_to_lambdify, fdx, 'numpy')


def gradient(func, x1, x2, point):
    fdx1 = fd(func, x1, (x1, x2))
    fdx2 = fd(func, x2, (x1, x2))
    return [fdx1(point[0], point[1]), fdx2(point[0], point[1])]


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]):
    betha0 = Symbol('betha0')
    f_func = lambdify((x1, x2), f, 'numpy')
    f_betha = f_func(x_point[0] + y_vector[0] * betha0, x_point[1] + y_vector[1] * betha0)
    f_betha = lambdify(betha0, f_betha, 'numpy')
    betha = 0
    betha_min = betha
    # Поиск минимума
    while betha <= 10:
        betha = betha + 1e-3
        if f_betha(betha_min) >= f_betha(betha):
            betha_min = betha
    # Значений для betha и g функций с ней
    betha_functions = []
    betha_variables = []
    i = 1
    for g in g_functions:
        # Переменная бета, т.е. бета1 бета2 бета3 и т.д.
        betha_variable = Symbol('betha' + str(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)
    result = np.array(result.args)
    min_result = min(result[0])
    if min_result < betha:
        return min_result
    return betha


def possible_directions_method(f: Function, g_functions, x0, epsilon):
    # Шаг 1
    x_ = np.array(x0)
    k = 0
    yk = []
    while k < 100:
        # Шаг 2
        f_gradient = gradient(f, x1, x2, x_)
        norm = 0
        for el in f_gradient:
            norm = norm + el**2
        if norm <= epsilon:

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

        else:
            # Шаг 3
            i_set = solve_system(g_functions, x_)
            # Шаг 4
            if len(i_set) == 0:
                yk = gradient(f, x1, x2, x_)
                yk = -1 * np.array(yk)
                zk = 2 * epsilon
            else:
                a = []
                b = [0]
                c = [0, 0, 1]  # значения коэффициентов функции для y1 y2 z
                row = gradient(f, x1, x2, x_)
                row.append(-1)  # задание значений для y1, y2 и -z
                a.append(row)
                for i in i_set:
                    if i != 0:
                        row = 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": False})
                # print('Результаты решения симплекс метода:\n', linprog_result)
                # шаг 5
                yk = [linprog_result.x[0], linprog_result.x[1]]
                zk = linprog_result.x[2]
            if yk[0] <= epsilon and yk[1] <= epsilon or m.fabs(zk) <= epsilon:

                # получение решения #
                return x_, yk, k, 5

            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
        print('Итерация', k)
    return x_, yk, k, -1


def main():
    x0 = [0, 0.95]
    epsilon = 0.03
    f = (x1 - 4) ** 2 + (x2 - 5) ** 2
    g = [x1 + x2 - 1, -x1, -x2]
    x_result, y_result, iter_num, exit_code = possible_directions_method(f, g, x0, epsilon)
    print("Задача была решена, точка x = ", x_result)
    if y_result is not None:
        print("Вектор y = ", y_result)
    print('Количество итераций: ', iter_num)
    print('Вышел после шага', exit_code)


if __name__ == "__main__":
    main()