Mở đầu

Trong bài viết dưới đây, chúng ta sẽ tìm hiểu cách sử dụng scipy để giải quyết những bài toán yêu cầu tối ưu hóa giá trị của một công thức toán học. Bài viết này sẽ không liên quan đến machine learning. Nhưng phương pháp đề cập trong bài viết này sẽ có một số ứng dụng thực tế. Bài viết này phù hợp với những người biết lập trình về python và phù hợp trong tất cả lĩnh vực.

Bài toán tối thiểu hóa

Chúng ta bắt đầu bằng bài toán đầu tiên.

Một công ty dầu có 2 nhà máy lọc dầu.
Nhà máy 1 mất \$20,000/ngày để chạy. Mỗi ngày nhà máy sẽ sản xuất được 400 thùng dầu cao cấp, 300 thùng dầu trung cấp, 200 thùng dầu cấp thấp.
Nhà máy 2 hiện đại hơn, mất \$25,000/ngày để chạy. Mỗi ngày nhà máy sản xuất được 300 thùng dầu cao cấp, 400 thùng dầu trung cấp, 500 thùng dầu cấp thấp mỗi ngày.
Công ty nhận được đơn hàng yêu cầu 25,000 thùng dầu cao cấp, 27,000 thùng dầu trung cấp, 30,000 thùng dầu cấp thấp.
Hỏi bao nhiêu ngày mỗi nhà máy lọc dầu cần chạy để giảm chi phí chạy đến mức tối thiểu mà vẫn có thể lọc dầu được đáp ứng được theo đơn hàng?

Gọi số ngày nhà máy 1 chạy là x1
Gọi số ngày nhà máy 2 chạy là x2
Tổng 2 nhà máy phải lọc được dầu với khối lượng lớn hơn hoặc bằng đơn hàng. Ta có các điều kiện sau:

  1. 400x1 + 300x2 >= 25,000
  2. 300x1 + 400x2 >= 27,000
  3. 200x1 + 500x2 >= 30,000
    Và mục tiêu của chúng ta là chi phí đến mức tối thiểu.
    Objective: Minimize 20000x1 + 25000x2
    Ta có thể tạo một chương trình để giải quyết vấn đề sau như thế nào?

Viết điều kiện dưới dạng code.

400x1 + 300x2 >= 25,000

sẽ được viết dưới dạng dictionary sau, thể hiện cho điều kiện

{"type": "ineq", "fun": lambda x: 400x[0] + 300x[1] - 25000}

type:"ineq" (inequality condition) có nghĩa là phương trình bất bình vì có dấu ">="
Scipy đòi hỏi tất cả các điều kiện phải được chuyển đổi về dạng sau:

f(x) >= c -> f(x) - c >= 0

f(x) <= c -> c - f(x) >= 0

Trong trường hợp:
400x[0] + 300x[1] >= 25000 sẽ được chuyển thành 400x[0] + 300x[1] - 25000 >= 0

và bạn chỉ cần để
{"fun": 400x[0] + 300x[1] - 25000}

import matplotlib.pyplot as plt
import numpy as np

constraints = [
    {"type": "ineq", "fun": lambda x: 400*x[0] + 300*x[1] - 25000},
    {"type": "ineq", "fun": lambda x: 300*x[0] + 400*x[1] - 27000},
    {"type": "ineq", "fun": lambda x: 200*x[0] + 500*x[1] - 30000},
    {"type": "ineq", "fun": lambda x: x[0]},
    {"type": "ineq", "fun": lambda x: x[1]}
]

Tạo một function objective thể hiện phương trình bạn muốn tối ưu hóa
Minimize 20000x1 + 25000x2

def objective(x):
    return (20000*x[0]) + (25000*x[1])

Giá trị bất kì ban đầu.

x0 = [0, 0]

from scipy.optimize import minimize
sol = minimize(
    objective, 
    x0,
    method = "COBYLA",
    constraints = constraints
)

Số ngày chạy nhà máy lọc dầu

sol.x

array([25., 50.])

Ta tìm thấy chi phí lọc dầu tối thiểu là:

round(objective(sol.x))

1750000.0

Bài toán tối đa hóa

Ứng dụng trong kinh tế: Chọn cách quảng cáo
Một công ty có lựa chọn quảng cáo thông qua TV, đài, và báo
Mỗi lần quảng cáo trên TV sẽ mất \$2000, trên báo mất \$600, trên đài mất \$300

Cho mỗi quảng cáo, TV sẽ quảng cáo được cho 100,000 người, báo cho 40,000 người, đài cho 18,000 người.
Trên báo, công ty sẽ chỉ được quảng cáo tối đa 10 lần một tuần.

Để cân bằng các phương tiện quảng cáo, công ty sẽ không được quảng cáo trên 50% cho đài, và phải quảng cáo ít nhất 10% trên TV.

Tổng ngân sách một tuần cho quảng cáo là $18,200.
Hỏi mỗi phương tiện sẽ có bao nhiêu quảng cáo, để có được tối đa số người xem quảng cáo.

Số quảng cáo cho TV, đài, báo là x1, x2, x3.

Tổng 2 nhà máy phải lọc được dầu với khối lượng lớn hơn hoặc bằng đơn hàng. Ta có các điều kiện sau:

  1. x1 + x2 + x3 <= 10
  2. 2000x1 + 600x2 + 300x3 <= 18200
  3. x2/(x1 + x2 + x3) <= 0.5
  4. x1/(x1 + x2 + x3) >= 0.1

Và mục tiêu của chúng ta là tối đa hóa.
Maximize 100000x1 + 40000x2 + 18000x3

Ta có thể tạo một chương trình để giải quyết vấn đề sau như thế nào?

constraints = [
    {"type": "ineq", "fun": lambda x: 10 - (x[0] + x[1] + x[2])},
    {"type": "ineq", "fun": lambda x: 18200 - (2000*x[0] + 600*x[1] + 300*x[2])},
    {"type": "ineq", "fun": lambda x: 0.5 - (x[1]/(x[0] + x[1] + x[2]))},
    {"type": "ineq", "fun": lambda x: 0.1 - (x[0]/(x[0] + x[1] + x[2]))},
    {"type": "ineq", "fun": lambda x: x[0]},
    {"type": "ineq", "fun": lambda x: x[1]},
    {"type": "ineq", "fun": lambda x: x[2]}
]

Một điểm cần lưu ý là, thay vì tối thiểu hóa, ta đang tối đa hóa phương trình mục tiêu. Ta có thể thay đổi tối đa hóa thành tối thiểu hóa bằng cách biến công thức mục tiêu thành âm. Vì tối thiểu hóa số âm là tối đa hóa số dương

def objective(x):
    return -(100000*x[0] + 40000*x[1] + 18000*x[2])

x0 = [0, 0, 0]

from scipy.optimize import minimize
sol = minimize(
    objective, 
    x0,
    method = "COBYLA",
    constraints = constraints
)

/home/leminhlong/.local/lib/python3.6/site-packages/ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in double_scalars
  after removing the cwd from sys.path.
/home/leminhlong/.local/lib/python3.6/site-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in double_scalars
  """

Số quảng cáo cho mỗi phương tiện (TV, đài, báo)

np.round(sol.x)

array([1., 5., 4.])

Số người xem tối đa là

-round(objective(sol.x))

372000.0

Cách tốt hơn để tìm giá trị thấp nhất

Trong những vấn đề về tối ưu, tìm giá trị ban đầu sẽ ảnh hưởng đến giá trị thấp nhất có thể ta tìm được. Global minima là giá trị nhỏ nhất thực tế có thể. Local minima là giá trị nhỏ nhất mà hiện tại ta tìm thấy được. Về lý tưởng chúng ta muốn tìm global minima.

Ghi chú: code được lấy từ https://www.youtube.com/watch?v=ms3aKKW_iRc

Tạo data để tối ưu hóa

import numpy as np
x = np.linspace(-10, 10, 51)
x

Phương trình để tối ưu hóa


```python
def f(x):
    return x**2 + 10*np.sin(x)

y = f(x)
y

fig = plt.figure(figsize=(5, 4), dpi=80)
plt.plot(x, y)
plt.grid('on')
plt.show()

/home/leminhlong/.local/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Passing one of 'on', 'true', 'off', 'false' as a boolean is deprecated; use an actual boolean (True/False) instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)

uc?id=1OgGaDtdmSvcwbWOxoMfvAAk7C7bRMrlz&export=download

x0 =8

Ta tìm giá trị tối thiểu của công thức, giống như hai bài toán trên

from scipy.optimize import minimize
sol = minimize(
    f, 
    x0=x0,
    method = "COBYLA",
)

float(sol.x)

3.8376023437500004

f(sol.x)

8.315585656391544

sol.x**2

14.727191748755496

Nếu sử dụng phương pháp tối ưu hóa thông thường như ở trên, từ điểm xuất phát ban đầu có giá trị là 78 ta đạt được giá trị thấp nhất là 8.3. Nhưng nếu nhìn kỹ trong hình có điểm thấp nhất nhỏ hơn 0, và code của chúng ta bị mắt kẹt tại đây. Điểm ta tìm được là local minima. Local minima ta tìm được phụ thuộc nhiều vào điểm ta bắt đầu

fig = plt.figure(figsize=(5, 4), dpi=80)
plt.plot(x, y)
plt.plot([x0, float(sol.x)], [f(x0), f(sol.x)], 'o')
plt.grid('on')
plt.show()

/home/leminhlong/.local/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Passing one of 'on', 'true', 'off', 'false' as a boolean is deprecated; use an actual boolean (True/False) instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)

uc?id=1NQv5ThvMOVRHc3EoHIR9fRchesk5fPKL&export=download

Có một phương pháp tốt hơn, bằng cách sử dụng basinhopping. Basinhopping sẽ cố gắng tìm giá trị thấp nhất của một phương trình, bằng cách tìm kiếm các điểm dựa vào các xác suất ngẫu nhiên lựa chọn liên tiếp (stochastic). Tuy vẫn có khả năng bashinhopping không tìm được điểm thấp nhất, nó vẫn có thể tốt hơn phương pháp tối thiểu hóa thông thường của ta ở trên.

T: temperature - giá trị càng cao thì độ nhảy từ điểm này đến điểm khác càng cao
stepsize - giá trị quan trọng trong bashinhopping. Một bước đi sẽ có giá trị (x0-stepsize) và (x0+stepsize) theo từng chiều của x0.

from scipy.optimize import basinhopping
sol = basinhopping(f, 
                 x0=6, 
                 T=1,  
                 stepsize=2)
sol

                        fun: -7.945823375615284
 lowest_optimization_result:       fun: -7.945823375615284
 hess_inv: array([[0.0858273]])
      jac: array([1.1920929e-07])
  message: 'Optimization terminated successfully.'
     nfev: 18
      nit: 5
     njev: 6
   status: 0
  success: True
        x: array([-1.30644001])
                    message: ['requested number of basinhopping iterations completed successfully']
      minimization_failures: 0
                       nfev: 1896
                        nit: 100
                       njev: 632
                          x: array([-1.30644001])

Sử dụng basinhopping. Ta đã tìm được điểm có giá trị tối thiểu của phương trình, global minima. Đối với các phương trình phức tạp hơn, có thể basinhopping vẫn chỉ tìm được local minima thôi

fig = plt.figure(figsize=(5, 4), dpi=80)
plt.plot(x, y)
plt.plot([x0, float(sol.x)], [f(x0), f(sol.x)], 'o')
plt.grid('on')
plt.show()

/home/leminhlong/.local/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Passing one of 'on', 'true', 'off', 'false' as a boolean is deprecated; use an actual boolean (True/False) instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)

uc?id=1D3Jvj3RuwclnhEeJVfzRQS4XXFBYwQNC&export=download

Kết luận

Tuy có giới hạn về các bài toán có thể áp dụng được. Scipy là một công cụ tốt cho những vấn đề tối ưu hóa công thức toán. Gần đây tôi có sử dụng phương pháp này để tối ưu hóa parameters cho công thức để gợi ý bài viết, mà không sử dụng machine learning. Hi vọng bài viết này sẽ giúp được trong công việc của các bạn trong lĩnh vực khác