Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import pylab
- from mpl_toolkits.mplot3d import Axes3D
- import numpy as np
- from scipy.optimize import minimize, fsolve, Bounds
- import matplotlib.pyplot as plt
- def makeData(): # данные для 3D графика
- x = np.arange(-5, 5, 0.1)
- y = np.arange(-5, 5, 0.1)
- xgrid, ygrid = np.meshgrid(x, y)
- zgrid = 3 * (ygrid + xgrid * xgrid)**2 + (ygrid * ygrid - 1)**2
- return xgrid, ygrid, zgrid
- def f(x): # целевая функция
- return 3 * (x[1] + x[0] * x[0])**2 + (x[1] * x[1] - 1)**2
- def system(x): # система для нахождения лямбды и точки минимума
- r = [0., 0., 0.]
- r[0] = 12*(x[1]+x[0]**2)*x[0]+4*(x[0]**2-1)*x[0]-0.5*x[2]
- r[1] = 6*(x[1]+x[0]**2) - x[2]
- r[2] = -(x[1]+0.5*x[0]+0.5)
- return r
- # отрисовка допустимой области
- figure, ax = plt.subplots()
- x1 = np.arange(-5, 5, 0.01)
- y1 = -0.5 * x1 - 0.5
- y3 = -0.8 + 0 * x1
- t = np.linspace(0, 2*np.pi, 100)
- plt.plot( -1 + np.sqrt(1.2) * np.cos(t) , 2 + np.sqrt(12) * np.sin(t) )
- plt.grid(color='lightgray',linestyle='--')
- ax.plot(x1, y1, linestyle='solid')
- ax.plot(x1, y3, linestyle='solid')
- for i in np.arange(-10,10,0.1):
- for j in np.arange(-10,10,0.1):
- if (j + 0.5*i+0.5) > 0 and j + 0.8 > 0 and -10*(i+1)**2-(j-2)**2+12 > 0:
- ax.scatter(i,j, s = 10, c = 'red')
- # отрисовка 3D графика
- fig = pylab.figure()
- axes = Axes3D(fig)
- x2, y2, z2 = makeData()
- axes.plot_surface(x2, y2, z2, alpha = 0.3)
- # применение minimize
- bnds = Bounds ([-np.inf, -np.inf], [np.inf, np.inf])
- ineq_cons = {'type': 'ineq',
- 'fun': lambda x: np.array ([(x[1]+0.5*x[0]+0.5),
- (-10*(x[0]+1)**2-(x[1]-2)**2+12),
- (x[1]+0.8)]),
- 'jac': lambda x: np.array ([[(0.5),(1.0)],
- [(-20*(x[0]+1)),(-2*(x[1]-2))],
- [(0.0),(1.0)]])
- }
- x0 = np.array([10.0, 10.0])
- res = minimize(f, x0, method='SLSQP',
- constraints=[ineq_cons], options={'ftol': 1e-9, 'disp': True},
- bounds=bnds)
- # рисуем получившуюся точку минимума на 3D графике
- axes.scatter(res.x[0], res.x[1], f(res.x), s = 5, c = 'red')
- # выводим справку-результат минимизации
- print('x = ' + str(res.x[0]) + ', y = ' + str(res.x[1]) + ', z(x,y) =' + str(f(res.x)))
- ax.scatter(res.x[0], res.x[1])
- # решаем систему (стационарность функции Лагранжа) и выводим результат
- ress = fsolve(system, (-1, 1, 1))
- print()
- print('After solution the system:')
- print(ress)
- # вывод графиков
- pylab.show()
- plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement