Advertisement
edward4324

Lab5

May 13th, 2022
641
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # численные методы поиска условного эсктремума
  2. # стр. 105
  3. import numpy as np
  4. import math as m
  5.  
  6. import scipy.optimize
  7. from sympy import *
  8. from scipy.optimize import linprog
  9.  
  10.  
  11. x1 = Symbol('x1')
  12. x2 = Symbol('x2')
  13. y1 = Symbol('y1')
  14. y2 = Symbol('y2')
  15. z = Symbol('z')
  16.  
  17.  
  18. def fd(func, x_to_diff, x_to_lambdify):
  19.     fdx = func.diff(x_to_diff)
  20.     return lambdify(x_to_lambdify, fdx, 'numpy')
  21.  
  22.  
  23. def gradient(func, x1, x2, point):
  24.     fdx1 = fd(func, x1, (x1, x2))
  25.     fdx2 = fd(func, x2, (x1, x2))
  26.     return [fdx1(point[0], point[1]), fdx2(point[0], point[1])]
  27.  
  28.  
  29. def solve_system(g_functions, x_point):
  30.     result = []
  31.     for i in range(0, len(g_functions)):
  32.         g = lambdify((x1, x2), g_functions[i], 'numpy')
  33.         if g(x_point[0], x_point[1]) == 0:
  34.             result.append(i)
  35.     return result
  36.  
  37.  
  38. def get_alpha_value(g_functions: list[Function], f: Function, x_point: list[float], y_vector: list[float]):
  39.     betha0 = Symbol('betha0')
  40.     f_func = lambdify((x1, x2), f, 'numpy')
  41.     f_betha = f_func(x_point[0] + y_vector[0] * betha0, x_point[1] + y_vector[1] * betha0)
  42.     f_betha = lambdify(betha0, f_betha, 'numpy')
  43.     betha = 0
  44.     betha_min = betha
  45.     # Поиск минимума
  46.     while betha <= 10:
  47.         betha = betha + 1e-3
  48.         if f_betha(betha_min) >= f_betha(betha):
  49.             betha_min = betha
  50.     # Значений для betha и g функций с ней
  51.     betha_functions = []
  52.     betha_variables = []
  53.     i = 1
  54.     for g in g_functions:
  55.         # Переменная бета, т.е. бета1 бета2 бета3 и т.д.
  56.         betha_variable = Symbol('betha' + str(i))
  57.         g_func = lambdify((x1, x2), g, 'numpy')
  58.         g_func = g_func(x_point[0] + y_vector[0] * betha_variable, x_point[1] + y_vector[1] * betha_variable)
  59.         betha_functions.append(g_func)
  60.         betha_variables.append(betha_variable)
  61.         i = i + 1
  62.     result = nonlinsolve(betha_functions, betha_variables)
  63.     result = np.array(result.args)
  64.     min_result = min(result[0])
  65.     if min_result < betha:
  66.         return min_result
  67.     return betha
  68.  
  69.  
  70. def possible_directions_method(f: Function, g_functions, x0, epsilon):
  71.     # Шаг 1
  72.     x_ = np.array(x0)
  73.     k = 0
  74.     yk = []
  75.     while k < 100:
  76.         # Шаг 2
  77.         f_gradient = gradient(f, x1, x2, x_)
  78.         norm = 0
  79.         for el in f_gradient:
  80.             norm = norm + el**2
  81.         if norm <= epsilon:
  82.  
  83.             # получение решения #
  84.             return x_, None, k, 2
  85.  
  86.         else:
  87.             # Шаг 3
  88.             i_set = solve_system(g_functions, x_)
  89.             # Шаг 4
  90.             if len(i_set) == 0:
  91.                 yk = gradient(f, x1, x2, x_)
  92.                 yk = -1 * np.array(yk)
  93.                 zk = 2 * epsilon
  94.             else:
  95.                 a = []
  96.                 b = [0]
  97.                 c = [0, 0, 1]  # значения коэффициентов функции для y1 y2 z
  98.                 row = gradient(f, x1, x2, x_)
  99.                 row.append(-1)  # задание значений для y1, y2 и -z
  100.                 a.append(row)
  101.                 for i in i_set:
  102.                     if i != 0:
  103.                         row = gradient(g_functions[i], x1, x2, x_)  # значения для y1 y2
  104.                         row.append(-1)  # значения для -z
  105.                         a.append(row)  # формирования коэффициентов для неравенств с j-ми функциями
  106.                         b.append(0)  # формирование ограничений для j-ых функций
  107.                 y1_bounds = [-1, 1]  # ограничения для y1
  108.                 y2_bounds = [-1, 1]  # ограничения для y2
  109.                 z_bounds = [None, None]
  110.                 bounds = [y1_bounds, y2_bounds, z_bounds]
  111.                 linprog_result = linprog(c, A_ub=a, b_ub=b, bounds=bounds, method='simplex', options={"disp": False})
  112.                 # print('Результаты решения симплекс метода:\n', linprog_result)
  113.                 # шаг 5
  114.                 yk = [linprog_result.x[0], linprog_result.x[1]]
  115.                 zk = linprog_result.x[2]
  116.             if yk[0] <= epsilon and yk[1] <= epsilon or m.fabs(zk) <= epsilon:
  117.  
  118.                 # получение решения #
  119.                 return x_, yk, k, 5
  120.  
  121.             else:
  122.                 # yk = [linprog_result.x[0], linprog_result.x[1]]
  123.                 # получение минимума из всех бета, что = альфа #
  124.                 alpha_k = get_alpha_value(g_functions, f, x_, yk)
  125.  
  126.                 # шаг 6
  127.                 # итерация #
  128.                 x_ = x_ + alpha_k * np.array(yk)
  129.         k = k + 1
  130.         print('Итерация', k)
  131.     return x_, yk, k, -1
  132.  
  133.  
  134. def main():
  135.     x0 = [0, 0.95]
  136.     epsilon = 0.03
  137.     f = (x1 - 4) ** 2 + (x2 - 5) ** 2
  138.     g = [x1 + x2 - 1, -x1, -x2]
  139.     x_result, y_result, iter_num, exit_code = possible_directions_method(f, g, x0, epsilon)
  140.     print("Задача была решена, точка x = ", x_result)
  141.     if y_result is not None:
  142.         print("Вектор y = ", y_result)
  143.     print('Количество итераций: ', iter_num)
  144.     print('Вышел после шага', exit_code)
  145.  
  146.  
  147. if __name__ == "__main__":
  148.     main()
  149.  
Advertisement
Advertisement
Advertisement
RAW Paste Data Copied
Advertisement