參考官網(wǎng)：Scipy.
對(duì)于有約束的最小化問題，Scipy提供的minimize這個(gè)包有三個(gè)：trust-constr, SLSQP'和COBYLA。它們要求使用稍微不同的結(jié)構(gòu)來定義約束。trust-constr需要要求約束被定義成一系列的LinearConstraint 和 NonlinearConstraint兩種類型。SLSQP'和COBYLA需要要求約束條件被定義為一連串的字典，其鍵為type、fun和jac。
考慮有約束的最小化2個(gè)變量的Rosenbrock函數(shù)
這個(gè)問題有唯一解：
[
x
0
,
x
1
]
=
[
0.4149
,
0.1701
]
[x_0,x_1]=[0.4149,0.1701]
[x0?,x1?]=[0.4149,0.1701]。
信賴域約束算法(method=‘trust-constr’)
信任域約束方法處理的是約束性最小化問題，其形式為:
除此之外，單邊約束可以通過對(duì)np.inf設(shè)置上限或下限并加上適當(dāng)?shù)姆?hào)來指定。
定義邊界約束
邊界限制 0 ≤ x 0 ≤ 1 , ? 0.5 ≤ x 1 ≤ 2.0 0≤x_0≤1,-0.5≤x_1≤2.0 0≤x0?≤1,?0.5≤x1?≤2.0被定義如下：
from scipy.optimize import Bounds
bounds = Bounds([0, -0.5], [1.0, 2.0])
定義線性約束
約束
x
0
+
2
x
1
≤
1
，
2
x
0
+
x
1
=
1
x_0+2x_1≤1，2x_0+x_1=1
x0?+2x1?≤1，2x0?+x1?=1可以寫成如下形式：
用LinearConstraint去定義：
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])
定義非線性約束
非線性約束為：
其雅可比矩陣(對(duì)每個(gè)變量求導(dǎo))為：
黑塞矩陣的線性組合：
用NonlinearConstraint去定義：
#定義非線性約束
def cons_f(x):
return [x[0]**2 + x[1], x[0]**2 - x[1]]
#定義導(dǎo)數(shù)
def cons_J(x):
return [[2*x[0], 1], [2*x[0], -1]]
#定義二階導(dǎo)數(shù)
def cons_H(x, v):
return v[0]*np.array([[2, 0], [0, 0]]) + v[1]*np.array([[2, 0], [0, 0]])
from scipy.optimize import NonlinearConstraint
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H)
另外，也可以將Hessian定義為一個(gè)稀疏矩陣。
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)
當(dāng)黑塞矩陣難以計(jì)算的時(shí)候，可以使用HessianUpdateStrategy，目前可用的策略有：BFGS和SR1。
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())
另外，Hessian可以用有限差分法進(jìn)行近似。
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess='2-point')
雅可比矩陣也可以用有限差分法估計(jì)，然而，在這種情況下，黑塞不能用有限差分計(jì)算，需要由用戶提供或用之前介紹的HessianUpdateStrategy定義：
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac='2-point', hess=BFGS())
求解
#設(shè)置初值
x0 = np.array([0.5, 0])
#這個(gè)需要結(jié)合之前無約束的算法來計(jì)算
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: 0.00e+00, execution time: 0.014 s.
[0.41494531 0.17010937]
完整代碼
import numpy as np
from scipy.optimize import minimize
#無約束
def rosen(x):
"""The Rosenbrock function"""
return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)
def rosen_der(x):
xm = x[1:-1]
xm_m1 = x[:-2]
xm_p1 = x[2:]
der = np.zeros_like(x)
der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm)
der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0])
der[-1] = 200*(x[-1]-x[-2]**2)
return der
def rosen_hess(x):
x = np.asarray(x)
H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
diagonal = np.zeros_like(x)
diagonal[0] = 1200*x[0]**2-400*x[1]+2
diagonal[-1] = 200
diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
H = H + np.diag(diagonal)
return H
#設(shè)置約束
from scipy.optimize import Bounds
bounds = Bounds([0, -0.5], [1.0, 2.0])
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])
def cons_f(x):
return [x[0]**2 + x[1], x[0]**2 - x[1]]
#定義導(dǎo)數(shù)
def cons_J(x):
return [[2*x[0], 1], [2*x[0], -1]]
#定義二階導(dǎo)數(shù)
def cons_H(x, v):
return v[0]*np.array([[2, 0], [0, 0]]) + v[1]*np.array([[2, 0], [0, 0]])
from scipy.optimize import NonlinearConstraint
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H)
#求解
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)
另外也可以對(duì)目標(biāo)函數(shù)的第一和第二導(dǎo)數(shù)進(jìn)行近似，例如，黑塞矩陣可以用SR1準(zhǔ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)
`gtol` termination condition is satisfied.
Number of iterations: 12, function evaluations: 24, CG iterations: 7, optimality: 4.48e-09, constraint violation: 0.00e+00, execution time: 0.02 s.
[0.41494531 0.17010937]
序列最小二乘(SLSQP) (method=‘SLSQP’)
SLSQP需要如下形式：
設(shè)置約束
線性和非線性約束都被定義為字典，其鍵為type、fun和jac：
首先是定義域：
from scipy.optimize import Bounds
bounds = Bounds([0, -0.5], [1.0, 2.0])
等式與不等式約束：文章來源：http://www.zghlxwxcb.cn/news/detail-788379.html
#不等式約束
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是對(duì)fun的求導(dǎo)
'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])}
求解
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 ]
完整代碼
import numpy as np
from scipy.optimize import minimize
#定義無約束的目標(biāo)函數(shù)
def rosen(x):
"""The Rosenbrock function"""
return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)
def rosen_der(x):
xm = x[1:-1]
xm_m1 = x[:-2]
xm_p1 = x[2:]
der = np.zeros_like(x)
der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm)
der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0])
der[-1] = 200*(x[-1]-x[-2]**2)
return der
def rosen_hess(x):
x = np.asarray(x)
H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
diagonal = np.zeros_like(x)
diagonal[0] = 1200*x[0]**2-400*x[1]+2
diagonal[-1] = 200
diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
H = H + np.diag(diagonal)
return H
#定義約束
from scipy.optimize import Bounds
bounds = Bounds([0, -0.5], [1.0, 2.0])
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])}
#求解
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)
PS：大部分 trust-constr方法可用的選項(xiàng)對(duì) SLSQP來說是不可用的。文章來源地址http://www.zghlxwxcb.cn/news/detail-788379.html
到了這里，關(guān)于【Scipy優(yōu)化使用教程】二、Scipy中有約束優(yōu)化的兩種算法的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
