03. 非线性规划
定义
如果目标函数或约束条件中包含非线性函数,就称这种规划问题为非线性规划问题。
非线性规划目前还没有适于各种问题的一般算法。
非线性规划模型描述如:
min
f
(
x
)
,
\min f(x),
minf(x),
s . t . { g i ( x ) ≤ 0 , i = 1 , 2 , … , m , h j ( x ) = 0 , j = 1 , 2 , … , l , s.t.\left\{\begin{aligned} &g_i(x)\le0,i=1,2,\dots,m,\\ &h_j(x)=0,j=1,2,\dots,l, \end{aligned}\right. s.t.{gi(x)≤0,i=1,2,…,m,hj(x)=0,j=1,2,…,l,
其中 x = [ x 1 , x 2 , … , x n ] T x=[x_1,x_2,\dots,x_n]^T x=[x1,x2,…,xn]T,而 f , g i , h j f,g_i,h_j f,gi,hj 都是实值函数。
一般非线性规划只能得到局部最优解,不能保证是全局最优解。
无约束非线性规划求解
设
f
(
x
)
f(x)
f(x) 具有连续的一阶偏导数,且
x
∗
x^*
x∗ 是无约束问题的局部极小点,则
∇
f
(
x
∗
)
=
0
\nabla f(x^*)=0
∇f(x∗)=0
其中
∇
f
(
x
)
\nabla f(x)
∇f(x) 表示
f
(
x
)
f(x)
f(x) 的梯度。
设
f
(
x
)
f(x)
f(x) 具有对各个变量的二阶偏导数,称矩阵
(
∂
2
f
∂
x
1
2
∂
2
f
∂
x
1
∂
x
2
…
∂
2
f
∂
x
1
∂
x
n
∂
2
f
∂
x
2
∂
x
1
∂
2
f
∂
x
2
2
…
∂
2
f
∂
x
2
∂
x
n
⋮
⋮
⋱
⋮
∂
2
f
∂
x
n
∂
x
1
∂
2
f
∂
x
n
∂
x
2
…
∂
2
f
∂
x
n
2
)
\begin{pmatrix} \frac{\partial^2f}{\partial x_1^2} & \frac{\partial^2f}{\partial x_1\partial x_2} & \dots & \frac{\partial^2f}{\partial x_1\partial x_n} \\ \frac{\partial^2f}{\partial x_2\partial x_1} & \frac{\partial^2f}{\partial x_2^2} & \dots & \frac{\partial^2f}{\partial x_2\partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2f}{\partial x_n\partial x_1} & \frac{\partial^2f}{\partial x_n\partial x_2} & \dots & \frac{\partial^2f}{\partial x_n^2} \\ \end{pmatrix}
⎝⎜⎜⎜⎜⎜⎛∂x12∂2f∂x2∂x1∂2f⋮∂xn∂x1∂2f∂x1∂x2∂2f∂x22∂2f⋮∂xn∂x2∂2f……⋱…∂x1∂xn∂2f∂x2∂xn∂2f⋮∂xn2∂2f⎠⎟⎟⎟⎟⎟⎞
为函数的黑塞矩阵,记
∇
2
f
(
x
)
\nabla^2f(x)
∇2f(x)。
因此,只要找到 x ∗ x^* x∗,满足 ∇ f ( x ∗ ) = 0 \nabla f(x^*)=0 ∇f(x∗)=0,且 ∇ 2 f ( x ∗ ) \nabla^2f(x^*) ∇2f(x∗) 为正定矩阵,则 x ∗ x^* x∗ 为无约束优化问题的局部最优解。
找到 x ∗ x^* x∗ 的具体方法有最速降线法、牛顿法等。
有约束非线性规划求解
有等式约束非线性规划的 Lagrange 乘数法
对于只有等式约束的非线性规划问题:
min
f
(
x
)
,
\min f(x),
minf(x),
s . t . { h j ( x ) = 0 , j = 1 , 2 , … , l , x ∈ R n . s.t.\left\{\begin{aligned} & h_j(x)=0,j=1,2,\dots,l,\\ & x\in R^n. \end{aligned}\right. s.t.{hj(x)=0,j=1,2,…,l,x∈Rn.
设函数
f
,
h
j
(
j
=
1
,
2
,
…
,
l
)
f,h_j(j=1,2,\dots,l)
f,hj(j=1,2,…,l) 在可行点
x
∗
x^*
x∗ 的某个邻域
N
(
x
∗
,
ϵ
)
N(x^*,\epsilon)
N(x∗,ϵ) 内可微,向量组
∇
h
j
(
x
∗
)
\nabla h_j(x^*)
∇hj(x∗) 线性无关,令
L
(
x
,
λ
)
=
f
(
x
)
−
λ
T
H
(
x
)
,
L(x,\lambda)=f(x)-\lambda^TH(x),
L(x,λ)=f(x)−λTH(x),
其中
λ
=
[
λ
1
,
λ
2
,
…
,
λ
l
]
T
\lambda=[\lambda_1,\lambda_2,\dots,\lambda_l]^T
λ=[λ1,λ2,…,λl]T
H ( x ) = [ h 1 ( x ) , h 2 ( x ) , … , h l ( x ) ] T H(x)=[h_1(x),h_2(x),\dots,h_l(x)]^T H(x)=[h1(x),h2(x),…,hl(x)]T
若
x
∗
x^*
x∗ 是局部最优解,则有
λ
∗
=
[
λ
1
∗
,
λ
2
∗
,
…
,
λ
l
∗
]
T
\lambda^*=[\lambda_1^*,\lambda_2^*,\dots,\lambda_l^*]^T
λ∗=[λ1∗,λ2∗,…,λl∗]T
使得
∇
L
(
x
∗
,
λ
∗
)
=
0
\nabla L(x^*,\lambda^*)=0
∇L(x∗,λ∗)=0,即
∇
f
(
x
∗
)
−
∑
j
=
1
l
λ
j
∗
∇
h
j
(
x
∗
)
=
0
\nabla f(x^*)-\sum_{j=1}^l\lambda_j^*\nabla h_j(x^*)=0
∇f(x∗)−j=1∑lλj∗∇hj(x∗)=0
从而将有约束条件转化为无约束条件。
有约束非线性规划的罚函数法
构造带参数的所谓增广目标函数,从而把有约束非线性规划问题转化为一系列无约束非线性规划问题。
增广目标函数通常由两部分组成:
- 原问题的目标函数;
- 约束条件构造出的惩罚项。
对于外点法,增广目标函数形式如下:
T
(
x
,
M
)
=
f
(
x
)
+
M
∑
i
=
1
m
[
max
{
0
,
g
i
(
x
)
}
]
+
M
∑
j
=
1
l
[
h
j
(
x
)
]
2
T(x,M)=f(x)+M\sum_{i=1}^m[\max\{0,g_i(x)\}]+M\sum_{j=1}^l[h_j(x)]^2
T(x,M)=f(x)+Mi=1∑m[max{0,gi(x)}]+Mj=1∑l[hj(x)]2
其中 M 是一个较大的正数。
从而转化为无约束问题:
min
T
(
x
,
M
)
,
x
∈
R
n
\min T(x,M),x\in R^n
minT(x,M),x∈Rn
罚函数法的计算精度可能较差。
Python 代码
利用 scipy、cvxopt、cvxpy 包,可以实现非线性规划求解。
scipy 求解
对于问题:
min
2
+
x
1
1
+
x
2
−
3
x
1
+
4
x
3
,
\min \frac{2+x_1}{1+x_2}-3x_1+4x_3,
min1+x22+x1−3x1+4x3,
s . t . 0.1 ≤ x i ≤ 0.9 , i = 1 , 2 , 3 s.t.\quad0.1\le x_i\le0.9,i=1,2,3 s.t.0.1≤xi≤0.9,i=1,2,3
使用 scipy 求解,代码如下:
# %%
import numpy as np
from scipy.optimize import minimize
# %%
def obj(x):
x1,x2,x3 = x
return (2+x1)/(1+x2)-3*x1+4*x3
bound = [(.1, .9) for _ in range(3)]
res = minimize(obj, np.ones(3), bounds=bound)
# %%
print('best obj =', res.fun)
print('success =', res.success)
print('best x =', res.x)
输出如下:
best obj = -0.7736842105263159
success = True
best x = [0.9 0.9 0.1]
cvxpy 求解
对于问题:
min
z
=
x
1
2
+
x
2
2
+
3
x
3
2
+
4
x
4
2
+
2
x
5
2
−
8
x
1
−
2
x
2
−
3
x
3
−
x
4
−
2
x
5
,
\min z=x_1^2+x_2^2+3x_3^2+4x_4^2+2x_5^2-8x_1-2x_2-3x_3-x_4-2x_5,
minz=x12+x22+3x32+4x42+2x52−8x1−2x2−3x3−x4−2x5,
s . t . { 0 ≤ x i ≤ 99 , i = 1 , 2 , … , 5 x i ∈ Z + , i = 1 , 2 , … , 5 x 1 + x 2 + x 3 + x 4 + x 5 ≤ 400 , x 1 + 2 x 2 + 2 x 3 + x 4 + 6 x 5 ≤ 800 , 2 x 1 + x 2 + 6 x 3 ≤ 200 , x 3 + x 4 + 5 x 5 ≤ 200. s.t.\left\{\begin{aligned} & 0\le x_i\le99,i=1,2,\dots,5\\ & x_i\in Z^+,i=1,2,\dots,5\\ & x_1+x_2+x_3+x_4+x_5\le400,\\ & x_1+2x_2+2x_3+x_4+6x_5\le800,\\ & 2x_1+x_2+6x_3\le200,\\ & x_3+x_4+5x_5\le200. \end{aligned}\right. s.t.⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧0≤xi≤99,i=1,2,…,5xi∈Z+,i=1,2,…,5x1+x2+x3+x4+x5≤400,x1+2x2+2x3+x4+6x5≤800,2x1+x2+6x3≤200,x3+x4+5x5≤200.
代码如下:
# %%
import numpy as np
import cvxpy as cp
# %%
# 目标函数的平方项
c1 = np.array([1, 1, 3, 4, 2])
# 目标函数的一次项
c2 = np.array([-8, -2, -3, -1, -2])
# 约束项
a = np.array([[1, 1, 1, 1, 1],
[1, 2, 2, 1, 6],
[2, 1, 6, 0, 0],
[0, 0, 1, 1, 5]])
b = np.array([400, 800, 200, 200])
# 决策变量
x = cp.Variable(5, integer=True)
# 目标函数
obj = cp.Minimize(c1 @ x**2 + c2 @ x)
# 约束
con = [0 <= x, x <= 99, a@x <= b]
# 问题模型
prob = cp.Problem(obj, con)
prob.solve(solver='CPLEX')
# %%
print('best z =', prob.value)
print('best x =', x.value)
输出如下:
best z = -17.0
best x = [4. 1. 0. 0. 0.]