作者：莫斑煒

編者按：本文使用SciPy的optimize模塊來求解非線性規劃問題，結合實際例子，引入非線性規劃問題的求解算法及相應函數的調用。

本文提綱

• 一維搜索／單變量優化問題

• 無約束多元優化問題

• 非線性最小二乘問題

• 約束優化問題

非線性規劃問題的目標函數或約束條件是非線性的。本文使用SciPy的optimize模塊來求解非線性規劃問題。

目標函數和約束條件是否連續光滑是非常重要的性質，這是因為如果光滑，則所有決策變量可微，多變量函數的偏導數組成的向量為梯度，梯度是指向目標函數增長最快的方向。將目標函數梯度作為搜索方向，對非線性規劃問題的求解具有重要的意義。這些函數或其導數\梯度的不連續性給許多現有的非線性優化問題的求解帶來了困難。在下文中，我們假設這些函數是連續且光滑的。

# Importing Modules from scipy import optimize import matplotlib.pyplot as plt import numpy as np import sympy

1、一維搜索／單變量優化問題（Univariate Optimization）

無約束非線性規劃最簡單的形式是一維搜索。一維搜索通常作為多維優化問題中的一部分出現，比如梯度下降法中每次最優迭代步長的估計。求解一維搜索常用的兩類方法是函數逼近法和區間收縮法。其中函數逼近法是指用較簡單的函數近似代替原來的函數，用近似函數的極小點來估計原函數的極小點，比如牛頓法；區間收縮法對于一個單谷函數通過迭代以不斷縮小該區間的長度，當區間長度足夠小時，可將該區間中的一點作為函數的極小點，比如黃金分割法。

e.g. 最小化一個單位體積的圓柱體的表面積。

Objective ? ? ? ? ? ? ?

? ?s.t. ? ? ? ? ? ? ? ? ? ? ?

目標函數為圓柱體的表面積，約束條件為圓柱體體積為1。變量為r、h，這是一個帶等式約束的二維優化問題。根據等式約束條件，可得 ? ? ? ? ? ? ?，代入目標函數得到一維搜索： ? ? ? ? ? ? ?。該問題可通過求解導數為0的方法來求解。

r, h = sympy.symbols("r, h") Area = 2 * sympy.pi * r**2 + 2 * sympy.pi * r * h Volume = sympy.pi * r**2 * h h_r = sympy.solve(Volume - 1)[0] Area_r = Area.subs(h_r) rsol = sympy.solve(Area_r.diff(r))[0] rsol.evalf() # 0.541926070139289# 再驗證二階導數為正且rsol對應最小值 Area_r.diff(r, 2).subs(r, rsol) Area_r.subs(r, rsol) _.evalf() # 5.53581044593209

使用SciPy中optimize模塊中的brent()函數來最小化一維函數。brent方法是一種混合方法，結合了牛頓法和二分法，既能保證穩定性又能快速收斂，通常是SciPy一維搜索的首選方法。

首先定義一個Python函數f作為目標函數。將f傳給optimize.brent，參數brack表示指定算法的開始區間。

def f(r):return 2 * np.pi * r**2 + 2 / r r_min = optimize.brent(f, brack=(0.1, 4)) r_min # 0.5419260772557135 f(r_min) # 5.535810445932086

通過以上求解結果可看出，r取值約為0.54，目標函數取值約為5.54。注意進行優化之前，最好通過目標函數可視化以確定合適的搜索區間或起點。

2、無約束多元優化問題

無約束多元優化問題的求解要比一維搜索問題的求解困難得多。多變量情況下，求解非線性梯度方程根的解析方法是不可行的，而黃金分割搜索法中使用的區間形式也不能直接用。最基本的方法是考慮目標函數f(x)在給定點x處的梯度。負梯度方向是指向函數f(x)減少最多的方向，即負梯度方向是函數的最速下降方向。牛頓多元優化方法是對最速下降法的一種改進，它可以提高收斂性。在單變量情況（一維搜索）下，牛頓法可以看作函數的局部二次逼近。在多變量情況下，與最速下降法相比，步長被函數Hessian矩陣的逆的所代替。

SciPy庫中，函數optimize.fmin_ncg使用牛頓法，optimize.fmin_ncg需要的參數包括目標函數和搜索起點，需要用于計算梯度的函數及用于計算Hessian矩陣的函數。

e.g. ? ? ? ? ? ? ?

使用牛頓法進行求解，需要知道梯度和Hessian矩陣，對每個變量使用sympy.diff函數（用于求導）以獲得梯度和Hessian矩陣。

x1, x2 = sympy.symbols("x_1, x_2") f_sym = (x1-1)**4 + 5 * (x2-1)**2- 2*x1*x2 fprime_sym = [f_sym.diff(x_) for x_ in (x1, x2)] # Gradient sympy.Matrix(fprime_sym) fhess_sym = [[f_sym.diff(x1_, x2_) for x1_ in (x1, x2)] for x2_ in (x1, x2)] # Hessian sympy.Matrix(fhess_sym)# 構造函數表達式 f_lmbda = sympy.lambdify((x1, x2), f_sym, 'numpy') fprime_lmbda = sympy.lambdify((x1, x2), fprime_sym, 'numpy') fhess_lmbda = sympy.lambdify((x1, x2), fhess_sym, 'numpy')def func_XY_to_X_Y(f):"""Wrapper for f(X) -> f(X[0], X[1])"""return lambda X: np.array(f(X[0], X[1])) f = func_XY_to_X_Y(f_lmbda) fprime = func_XY_to_X_Y(fprime_lmbda) fhess = func_XY_to_X_Y(fhess_lmbda) # f、fprime、fhess不僅作用于單個值，還作用于整個向量# Newton method x_opt = optimize.fmin_ncg(f, (0, 0), fprime=fprime, fhess=fhess) # (0,0)為搜索起點 x_opt # array([1.88292613, 1.37658523])

? ? ? ?

目標函數極小值用紅星標出

由于牛頓法需要提供梯度和Hessian矩陣的計算函數，有時它們很難計算。下面介紹的方法可以不用傳入計算梯度和Hessian矩陣的函數。用optimize.fmin_bfgs（BFGS算法）和optimize.fmin_cg（共軛梯度法）求解，不必為Hessian函數提供函數。如果沒有計算梯度的函數也可以不傳入，直接調用。

# 僅提供梯度，不提供hessian矩陣 # BFGS method x_opt = optimize.fmin_bfgs(f, (0, 0), fprime=fprime) x_opt# conjugate gradient method x_opt = optimize.fmin_cg(f, (0, 0), fprime=fprime) x_opt

也可以在不提供梯度函數的情況下使用：

x_opt = optimize.fmin_bfgs(f, (0, 0))

但是如果傳遞了計算梯度的函數，求解結果將更好，求解速度也更快。

在梯度和Hessian矩陣都未知的情況下，可使用BFGS算法（擬牛頓方法）。如果梯度和Hessian矩陣都是已知的，那么牛頓法是通常收斂速度最快的方法，但是對于大規模問題牛頓法的每次迭代的計算量較大。雖然BFGS算法和共軛梯度法在理論上的收斂速度比牛頓法慢，但它們有時穩定性更好。

以上介紹的多元優化方法只能收斂到局部最優，即使存在全局最優解，這些方法也很難跳出局部最優，對于這種情況可以通過蠻力搜索找到更好的搜索起點來提升算法的性能。SciPy提供optimize.brute函數執行這種系統搜索。

e.g. ? ? ? ? ? ? ?

def f(X):x, y = Xreturn (4 * np.sin(np.pi * x) + 6 * np.sin(np.pi * y)) + (x- 1)**2 + (y- 1)**2

調用optimize.brute，目標函數f作為第一個參數，切片對象元組作為第二個參數，每個坐標一個。切片對象指定了搜索范圍，并通過設置參數finish=none直接給出蠻力的搜索方案。

x_start = optimize.brute(f, (slice(-3, 5, 0.5), slice(-3, 5, 0.5)), finish=None) x_start # array([2. , 1.5]) f(x_start) # 得到更好的搜索起點后代入optimize.fmin_bfgs中 x_opt = optimize.fmin_bfgs(f, x_start) # 可以和任意搜索起點作對比 x_opt = optimize.fmin_bfgs(f,(0,0))

可以看出目標函數更優，且求解時間更短，迭代次數更少，速度更快。

minimize()函數為scipy.optimize中的多變量標量函數提供了無約束和約束最小化算法的通用接口。SciPy提供為多變量優化求解器提供統一的接口optimize.minimize()，它通過method參數為求解器指定特定的函數，這樣可以更容易地在不同的求解器之間切換。一維搜索的接口為optimize.scalar_minimize。

x_opt = optimize.fmin_bfgs(f, x_start) # 等價于下面的語句 result = optimize.minimize(f, x_start, method= 'BFGS') # 也可不指定method x_opt = result.x

3、非線性最小二乘問題

? ? ? ? ? ? ? ? ? ?

optimize.leatsq函數用于求解最佳最小二乘擬合。

e.g. ? ? ? ? ? ? ?

# 觀測值 beta = (0.25, 0.75, 0.5) def f(x, b0, b1, b2):return b0 + b1 * np.exp(-b2 * x**2) xdata = np.linspace(0, 5, 50) y = f(xdata, *beta) # 加入噪音 ydata = y + 0.05 * np.random.randn(len(xdata))def g(beta):return ydata- f(xdata, *beta) beta_start = (1, 1, 1) beta_opt, beta_cov = optimize.leastsq(g, beta_start) beta_opt # array([0.24935676, 0.74672532, 0.49918151]) # 求解結果接近beta(0.25, 0.75, 0.5)

擬合曲線：

? ? ? ? ? ? ?

SciPy的optimize模塊還通過函數optimize.curve_fit為非線性最小二乘擬合提供了另一種接口，它不用向函數傳遞一個帶剩余誤差(residuals)的顯式函數。

beta_opt, beta_cov = optimize.curve_fit(f, xdata, ydata) beta_opt # array([0.24935676, 0.74672532, 0.49918151])

4、約束優化問題

約束條件增加了非線性規劃問題求解的復雜性。可以使用拉格朗日乘子法（Lagrange multiplier），通過引入附加變量將約束優化問題轉化為無約束問題。

e.g.

? ? ? ? ? ? ?

求構成的最大體積。

求解該問題使用拉格朗日乘子法，拉格朗日函數為 ? ? ? ? ? ? ?。可以找到使得該函數梯度為0的駐點求解該問題。

x = x0, x1, x2, l = sympy.symbols("x_0, x_1, x_2, lambda") f = x0 * x1 * x2 g = 2 * (x0 * x1 + x1 * x2 + x2 * x0)- 1 L = f + l * g grad_L = [sympy.diff(L, x_) for x_ in x] sols = sympy.solve(grad_L) sols # 求解結果中負值舍去 g.subs(sols[0]) # 0 f.subs(sols[0])

SciPy的optimize模塊中的函數slsqp()或通過optimize.minimize()函數設置method為'SLSQP'，為了用SciPy的slsqp求解器求解問題，我們需要為目標函數和約束函數定義python函數。由于optimize模塊優化函數用于解決最小化問題，而這里是求解最大化問題，故將函數f設置為原始目標函數的相反數。

def f(X):return -X[0] * X[1] * X[2] def g(X):return 2 * (X[0]*X[1] + X[1] * X[2] + X[2] * X[0])- 1

接下來定義g(x)=0的約束。通過字典形式傳入約束條件，該字典中key（value）是type（'eq'或'ineq'）、fun（約束條件）、jac（約束函數的Jacobian矩陣）和args（約束函數的其他參數和計算其Jacobian的函數）。最后調用optimize.minimize函數。

constraint = dict(type='eq', fun=g) result = optimize.minimize(f, [0.5,1,1.5],method='SLSQP',constraints=[constraint]) result

可以看出通過以上方法求得的解與通過使用拉格朗日乘子法得到的解基本一致。

要解決不等式約束問題，只需在約束字典中設置type='ineq'，并給出相應的不等式函數。

e.g.

? ? ? ? ? ? ?

def f(X):return (X[0]- 1)**2 + (X[1]- 1)**2 def g(X):return X[1]- 1.75- (X[0]- 0.75)**4 constraints = [dict(type='ineq', fun=g)] x_opt = optimize.minimize(f, (0, 0), method='BFGS').x #無約束 x_cons_opt = optimize.minimize(f(0,0),method='SLSQP',constraints=constraints).x

? ? ???

約束問題的可行域為灰色陰影，紅星和藍星分別為有約束和無約束問題的最優解。

參考文獻：

Johansson R. Numerical Python: Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib[M]. Apress, 2019.

備注：公眾號菜單包含了整理了一本AI小抄，非常適合在通勤路上用學習。

往期精彩回顧2019年公眾號文章精選適合初學者入門人工智能的路線及資料下載機器學習在線手冊深度學習在線手冊AI基礎下載（第一部分）備注：加入本站微信群或者qq群，請回復“加群”加入知識星球（4500+用戶，ID：92416895），請回復“知識星球”

喜歡文章，點個在看

總結

以上是生活随笔為你收集整理的优化 | 利用SciPy求解非线性规划问题的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。