作者：莫斑煒

編者按：本文使用SciPy的optimize模塊來求解非線性規(guī)劃問題，結(jié)合實際例子，引入非線性規(guī)劃問題的求解算法及相應(yīng)函數(shù)的調(diào)用。

本文提綱

• 一維搜索／單變量優(yōu)化問題

• 無約束多元優(yōu)化問題

• 非線性最小二乘問題

• 約束優(yōu)化問題

非線性規(guī)劃問題的目標(biāo)函數(shù)或約束條件是非線性的。本文使用SciPy的optimize模塊來求解非線性規(guī)劃問題。

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

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

1、一維搜索／單變量優(yōu)化問題（Univariate Optimization）

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

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

Objective ? ? ? ? ? ? ?

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

目標(biāo)函數(shù)為圓柱體的表面積，約束條件為圓柱體體積為1。變量為r、h，這是一個帶等式約束的二維優(yōu)化問題。根據(jù)等式約束條件，可得 ? ? ? ? ? ? ?，代入目標(biāo)函數(shù)得到一維搜索： ? ? ? ? ? ? ?。該問題可通過求解導(dǎo)數(shù)為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# 再驗證二階導(dǎo)數(shù)為正且rsol對應(yīng)最小值 Area_r.diff(r, 2).subs(r, rsol) Area_r.subs(r, rsol) _.evalf() # 5.53581044593209

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

首先定義一個Python函數(shù)f作為目標(biāo)函數(shù)。將f傳給optimize.brent，參數(shù)brack表示指定算法的開始區(qū)間。

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

通過以上求解結(jié)果可看出，r取值約為0.54，目標(biāo)函數(shù)取值約為5.54。注意進(jìn)行優(yōu)化之前，最好通過目標(biāo)函數(shù)可視化以確定合適的搜索區(qū)間或起點。

2、無約束多元優(yōu)化問題

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

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

e.g. ? ? ? ? ? ? ?

使用牛頓法進(jìn)行求解，需要知道梯度和Hessian矩陣，對每個變量使用sympy.diff函數(shù)（用于求導(dǎo)）以獲得梯度和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)# 構(gòu)造函數(shù)表達(dá)式 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])

? ? ? ?

目標(biāo)函數(shù)極小值用紅星標(biāo)出

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

# 僅提供梯度，不提供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

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

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

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

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

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

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

調(diào)用optimize.brute，目標(biāo)函數(shù)f作為第一個參數(shù)，切片對象元組作為第二個參數(shù)，每個坐標(biāo)一個。切片對象指定了搜索范圍，并通過設(shè)置參數(shù)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))

可以看出目標(biāo)函數(shù)更優(yōu)，且求解時間更短，迭代次數(shù)更少，速度更快。

minimize()函數(shù)為scipy.optimize中的多變量標(biāo)量函數(shù)提供了無約束和約束最小化算法的通用接口。SciPy提供為多變量優(yōu)化求解器提供統(tǒng)一的接口optimize.minimize()，它通過method參數(shù)為求解器指定特定的函數(shù)，這樣可以更容易地在不同的求解器之間切換。一維搜索的接口為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函數(shù)用于求解最佳最小二乘擬合。

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]) # 求解結(jié)果接近beta(0.25, 0.75, 0.5)

擬合曲線：

? ? ? ? ? ? ?

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

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

4、約束優(yōu)化問題

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

e.g.

? ? ? ? ? ? ?

求構(gòu)成的最大體積。

求解該問題使用拉格朗日乘子法，拉格朗日函數(shù)為 ? ? ? ? ? ? ?。可以找到使得該函數(shù)梯度為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 # 求解結(jié)果中負(fù)值舍去 g.subs(sols[0]) # 0 f.subs(sols[0])

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

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（約束函數(shù)的Jacobian矩陣）和args（約束函數(shù)的其他參數(shù)和計算其Jacobian的函數(shù)）。最后調(diào)用optimize.minimize函數(shù)。

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

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

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

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

? ? ???

約束問題的可行域為灰色陰影，紅星和藍(lán)星分別為有約束和無約束問題的最優(yōu)解。

參考文獻(xiàn)：

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

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

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

喜歡文章，點個在看

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。