SciPy (phát âm là sai pie) là một gói toán học dựa trên NumPy cũng bao gồm các thư viện C và Fortran. SciPy biến phiên Python tương tác của bạn thành môi trường khoa học dữ liệu hoàn chỉnh như MATLAB, IDL, Octave, R hoặc SciLab.
Trong bài viết này, chúng ta sẽ xem xét các kỹ thuật cơ bản của lập trình toán học - giải các bài toán tối ưu hóa có điều kiện cho hàm vô hướng nhiều biến bằng cách sử dụng gói scipy.optimize. Các thuật toán tối ưu hóa không giới hạn đã được thảo luận trong . Bạn luôn có thể nhận được trợ giúp chi tiết và cập nhật hơn về các hàm scipy bằng cách sử dụng lệnh help(), Shift+Tab hoặc trong .
Giới thiệu
Một giao diện chung để giải quyết cả các vấn đề tối ưu hóa có điều kiện và không bị ràng buộc trong gói scipy.optizes được cung cấp bởi hàm minimize(). Tuy nhiên, người ta biết rằng không có một phương pháp chung nào để giải quyết mọi vấn đề, vì vậy việc lựa chọn một phương pháp phù hợp như mọi khi đều thuộc về người nghiên cứu.
Thuật toán tối ưu hóa thích hợp được chỉ định bằng cách sử dụng đối số hàm minimize(..., method="").
Để tối ưu hóa có điều kiện một hàm nhiều biến, có thể triển khai các phương pháp sau:
trust-constr- tìm kiếm mức tối thiểu cục bộ trong vùng tin cậy. , ;SLSQP— lập trình bậc hai tuần tự có ràng buộc, phương pháp Newton để giải hệ Lagrange. .TNC- Cắt ngắn ràng buộc Newton, hạn chế số lần lặp, tốt cho các hàm phi tuyến có số lượng biến độc lập lớn. .L-BFGS-B— một phương pháp của nhóm Broyden–Fletcher–Goldfarb–Shanno, được triển khai với mức tiêu thụ bộ nhớ giảm do tải một phần vectơ từ ma trận Hessian. , .COBYLA— Tối ưu hóa ràng buộc MARE bằng xấp xỉ tuyến tính, tối ưu hóa ràng buộc với xấp xỉ tuyến tính (không tính toán độ dốc). .
Tùy thuộc vào phương pháp đã chọn, các điều kiện và hạn chế để giải quyết vấn đề được đặt ra khác nhau:
- đối tượng lớp
Boundsđối với các phương thức L-BFGS-B, TNC, SLSQP, Trust-constr; - danh sách
(min, max)đối với các phương thức tương tự L-BFGS-B, TNC, SLSQP, Trust-constr; - một đối tượng hoặc một danh sách các đối tượng
LinearConstraint,NonlinearConstraintđối với các phương thức COBYLA, SLSQP, Trust-constr; - từ điển hoặc danh sách từ điển
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}cho các phương pháp COBYLA, SLSQP.
Đề cương bài viết:
1) Xem xét việc sử dụng thuật toán tối ưu hóa có điều kiện trong vùng tin cậy (method=”trust-constr”) với các ràng buộc được chỉ định làm đối tượng Bounds, LinearConstraint, NonlinearConstraint ;
2) Xem xét lập trình tuần tự bằng phương pháp bình phương nhỏ nhất (phương thức = "SLSQP") với các hạn chế được chỉ định dưới dạng từ điển {'type', 'fun', 'jac', 'args'};
3) Phân tích một ví dụ về tối ưu hóa các sản phẩm được sản xuất bằng cách sử dụng ví dụ về một studio trên web.
Phương pháp tối ưu hóa có điều kiện="trust-constr"
Thực hiện phương pháp trust-constr dựa trên cho các vấn đề có ràng buộc về hình thức bình đẳng và về cho các bài toán có ràng buộc ở dạng bất đẳng thức. Cả hai phương pháp đều được thực hiện bằng thuật toán để tìm cực tiểu cục bộ trong vùng tin cậy và rất phù hợp cho các bài toán quy mô lớn.
Công thức toán học của bài toán tìm cực tiểu ở dạng tổng quát:
Đối với các ràng buộc bình đẳng nghiêm ngặt, giới hạn dưới được đặt bằng giới hạn trên .
Đối với ràng buộc một chiều, giới hạn trên hoặc giới hạn dưới được đặt np.inf bằng dấu tương ứng.
Cần tìm giá trị nhỏ nhất của hàm Rosenbrock đã biết của hai biến:
Trong trường hợp này, các hạn chế sau được đặt trên miền định nghĩa của nó:
Trong trường hợp của chúng tôi, có một giải pháp duy nhất tại thời điểm , trong đó chỉ hạn chế thứ nhất và thứ tư là hợp lệ.
Chúng ta hãy xem xét các hạn chế từ dưới lên trên và xem cách chúng ta có thể viết chúng trong scipy.
Hạn chế и hãy xác định nó bằng đối tượng Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Hạn chế и Hãy viết nó ở dạng tuyến tính:
Hãy xác định những ràng buộc này như một đối tượng LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])Và cuối cùng là ràng buộc phi tuyến ở dạng ma trận:
Chúng tôi xác định ma trận Jacobian cho ràng buộc này và sự kết hợp tuyến tính của ma trận Hessian với một vectơ tùy ý :
Bây giờ chúng ta có thể định nghĩa một ràng buộc phi tuyến là một đối tượng NonlinearConstraint:
from scipy.optimize import NonlinearConstraint
def cons_f(x):
return [x[0]**2 + x[1], x[0]**2 - x[1]]
def cons_J(x):
return [[2*x[0], 1], [2*x[0], -1]]
def cons_H(x, v):
return v[0]*np.array([[2, 0], [0, 0]]) + v[1]*np.array([[2, 0], [0, 0]])
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H)Nếu kích thước lớn, ma trận cũng có thể được chỉ định ở dạng thưa thớt:
from scipy.sparse import csc_matrix
def cons_H_sparse(x, v):
return v[0]*csc_matrix([[2, 0], [0, 0]]) + v[1]*csc_matrix([[2, 0], [0, 0]])
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
jac=cons_J, hess=cons_H_sparse)hoặc như một đối tượng LinearOperator:
from scipy.sparse.linalg import LinearOperator
def cons_H_linear_operator(x, v):
def matvec(p):
return np.array([p[0]*2*(v[0]+v[1]), 0])
return LinearOperator((2, 2), matvec=matvec)
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
jac=cons_J, hess=cons_H_linear_operator)Khi tính ma trận Hessian đòi hỏi nhiều nỗ lực, bạn có thể sử dụng một lớp . Các chiến lược sau đây có sẵn: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hessian cũng có thể được tính bằng cách sử dụng sai phân hữu hạn:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Ma trận Jacobian cho các ràng buộc cũng có thể được tính bằng cách sử dụng sai phân hữu hạn. Tuy nhiên, trong trường hợp này, ma trận Hessian không thể tính được bằng sai phân hữu hạn. Hessian phải được định nghĩa là một hàm hoặc sử dụng lớp HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())Giải pháp cho vấn đề tối ưu hóa trông như thế này:
from scipy.optimize import minimize
from scipy.optimize import rosen, rosen_der, rosen_hess, rosen_hess_prod
x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess,
constraints=[linear_constraint, nonlinear_constraint],
options={'verbose': 1}, bounds=bounds)
print(res.x)`gtol` termination condition is satisfied.
Number of iterations: 12, function evaluations: 8, CG iterations: 7, optimality: 2.99e-09, constraint violation: 1.11e-16, execution time: 0.033 s.
[0.41494531 0.17010937]Nếu cần, hàm tính toán Hessian có thể được xác định bằng lớp LinearOperator
def rosen_hess_linop(x):
def matvec(p):
return rosen_hess_prod(x, p)
return LinearOperator((2, 2), matvec=matvec)
res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess_linop,
constraints=[linear_constraint, nonlinear_constraint],
options={'verbose': 1}, bounds=bounds)
print(res.x)hoặc tích của Hessian và một vectơ tùy ý thông qua tham số hessp:
res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hessp=rosen_hess_prod,
constraints=[linear_constraint, nonlinear_constraint],
options={'verbose': 1}, bounds=bounds)
print(res.x)Ngoài ra, đạo hàm bậc nhất và bậc hai của hàm đang được tối ưu hóa có thể được tính gần đúng. Ví dụ: Hessian có thể được tính gần đúng bằng cách sử dụng hàm SR1 (xấp xỉ gần đúng Newton). Độ dốc có thể được tính gần đúng bằng sai phân hữu hạn.
from scipy.optimize import SR1
res = minimize(rosen, x0, method='trust-constr', jac="2-point", hess=SR1(),
constraints=[linear_constraint, nonlinear_constraint],
options={'verbose': 1}, bounds=bounds)
print(res.x)Phương pháp tối ưu hóa có điều kiện="SLSQP"
Phương pháp SLSQP được thiết kế để giải quyết các bài toán cực tiểu hóa hàm có dạng:
Где и - tập hợp các chỉ số biểu thức mô tả các hạn chế dưới dạng đẳng thức hoặc bất đẳng thức. - tập hợp các giới hạn dưới và trên cho miền định nghĩa của hàm.
Các ràng buộc tuyến tính và phi tuyến được mô tả dưới dạng từ điển với các phím type, fun и jac.
ineq_cons = {'type': 'ineq',
'fun': lambda x: np.array ([1 - x [0] - 2 * x [1],
1 - x [0] ** 2 - x [1],
1 - x [0] ** 2 + x [1]]),
'jac': lambda x: np.array ([[- 1.0, -2.0],
[-2 * x [0], -1.0],
[-2 * x [0], 1.0]])
}
eq_cons = {'type': 'eq',
'fun': lambda x: np.array ([2 * x [0] + x [1] - 1]),
'jac': lambda x: np.array ([2.0, 1.0])
}Việc tìm kiếm giá trị nhỏ nhất được thực hiện như sau:
x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='SLSQP', jac=rosen_der,
constraints=[eq_cons, ineq_cons], options={'ftol': 1e-9, 'disp': True},
bounds=bounds)
print(res.x)Optimization terminated successfully. (Exit mode 0)
Current function value: 0.34271757499419825
Iterations: 4
Function evaluations: 5
Gradient evaluations: 4
[0.41494475 0.1701105 ]Ví dụ về tối ưu hóa
Liên quan đến việc chuyển đổi sang cơ cấu công nghệ thứ năm, chúng ta hãy xem xét việc tối ưu hóa sản xuất bằng ví dụ về một studio web, mang lại cho chúng ta một khoản thu nhập nhỏ nhưng ổn định. Hãy tưởng tượng mình là giám đốc của một nhà bếp sản xuất ba loại sản phẩm:
- x0 - bán landing page, từ 10 tr.
- x1 - website công ty, từ 20 tr.
- x2 - cửa hàng trực tuyến, từ 30 tr.
Đội ngũ làm việc thân thiện của chúng tôi bao gồm bốn cấp dưới, hai cấp trung và một cấp cao. Quỹ thời gian làm việc hàng tháng của họ:
- Tháng Sáu:
4 * 150 = 600 чел * час, - giữa:
2 * 150 = 300 чел * час, - thưa ông:
150 чел * час.
Hãy để cấp dưới sẵn có đầu tiên dành (0, 1, 2) giờ cho việc phát triển và triển khai một trang loại (x10, x20, x30), cấp trung - (7, 15, 20), cấp cao - (5, 10, 15 ) giờ trong khoảng thời gian đẹp nhất trong cuộc đời bạn.
Giống như bất kỳ giám đốc bình thường nào, chúng tôi muốn tối đa hóa lợi nhuận hàng tháng. Bước đầu tiên để thành công là viết ra hàm mục tiêu value là số tiền thu nhập từ sản phẩm được sản xuất mỗi tháng:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Đây không phải là lỗi; khi tìm giá trị cực đại, hàm mục tiêu được cực tiểu hóa với dấu ngược lại.
Bước tiếp theo là cấm nhân viên của chúng tôi làm việc quá sức và đưa ra các hạn chế về giờ làm việc:
Tương đương là gì:
ineq_cons = {'type': 'ineq',
'fun': lambda x: np.array ([600 - 10 * x [0] - 20 * x [1] - 30 * x[2],
300 - 7 * x [0] - 15 * x [1] - 20 * x[2],
150 - 5 * x [0] - 10 * x [1] - 15 * x[2]])
}Một hạn chế chính thức là sản lượng sản phẩm chỉ phải dương:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])Và cuối cùng, giả định màu hồng nhất là vì giá rẻ và chất lượng cao nên hàng dài khách hàng hài lòng liên tục xếp hàng chờ chúng tôi. Chúng ta có thể tự mình chọn khối lượng sản xuất hàng tháng, dựa trên việc giải quyết vấn đề tối ưu hóa bị ràng buộc với scipy.optimize:
x0 = np.array([10, 10, 10])
res = minimize(value, x0, method='SLSQP', constraints=ineq_cons, bounds=bnds)
print(res.x)[7.85714286 5.71428571 3.57142857]Hãy làm tròn số thành số nguyên và tính khối lượng hàng tháng của người chèo thuyền với mức phân bổ sản phẩm tối ưu x = (8, 6, 3) :
- Tháng Sáu:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - giữa:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - thưa ông:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Kết luận: để giám đốc nhận được mức tối đa xứng đáng của mình, tối ưu nhất là tạo 8 trang đích, 6 trang cỡ trung bình và 3 cửa hàng mỗi tháng. Trong trường hợp này, đàn anh phải cày mà không được nhìn lên khỏi máy, tải trọng của đàn giữa sẽ xấp xỉ 2/3, đàn em ít hơn một nửa.
Kết luận
Bài viết nêu ra những kỹ thuật cơ bản để làm việc với gói scipy.optimize, được sử dụng để giải các bài toán tối thiểu hóa có điều kiện. Cá nhân tôi sử dụng scipy hoàn toàn vì mục đích học thuật, đó là lý do tại sao ví dụ được đưa ra lại mang tính chất hài hước như vậy.
Bạn có thể tìm thấy rất nhiều lý thuyết và ví dụ ảo, chẳng hạn như trong cuốn sách “Lập trình toán học trong các ví dụ và bài toán” của I.L. Akulich. Nhiều ứng dụng khó tính hơn scipy.optimize để xây dựng cấu trúc 3D từ một tập hợp hình ảnh () có thể được xem trong .
Nguồn thông tin chính là những người muốn đóng góp vào việc dịch phần này và phần khác scipy Chào mừng bạn đến .
Cảm ơn tham gia chuẩn bị xuất bản.
Nguồn: www.habr.com