# 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
78.         norm = 0
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.