一般的非线性规划求解(非凸函数)
记录奋斗路上的点滴，
希望能帮到一样刻苦的你！
如有不足欢迎指正！
共同学习交流！
🌎欢迎各位→点赞 👍+ 收藏⭐ + 留言📝
每一个裂缝都是为透出光而努力！
一起加油！
一般的非线性规划求解(非凸函数)
例 5.5 (续例 5.3) 求解如下二次规划模型
min−x12−0.3x1x2−2x22+98x1+277x2s.t.{x1+x2≤100,x1−2x2≤0,x1,x2≥0.
\begin{aligned}
& \min -x_1^2 - 0.3 x_1 x_2 - 2 x_2^2 + 98 x_1 + 277 x_2 \\
& \text{s.t.} \begin{cases}
x_1 + x_2 \leq 100, \\
x_1 - 2 x_2 \leq 0, \\
x_1, x_2 \geq 0.
\end{cases}
\end{aligned}
min−x12−0.3x1x2−2x22+98x1+277x2s.t.⎩⎨⎧x1+x2≤100,x1−2x2≤0,x1,x2≥0.
由于上述二次规划的目标函数不是凸函数，不能使用 cvxpy 库进行求解。
使用 Python 软件求得的最优解为
x1=0,x2=0,
x_1 = 0, \quad x_2 = 0,
x1=0,x2=0,
目标函数的最优值为 0。
计算的 Python 程序如下：
from scipy.optimize import minimize
import numpy as np# 定义目标向量
c1 = np.array([[-1, -0.15], [-0.15, -2]])
c2 = np.array([98, 277])
# 创建目标函数
obj = lambda x: x @ c1 @ x + c2 @ x
d = np.array([[1, 1], [1, -2]])
e = np.array([100, 0])
# 创建约束条件
constraints = {'type': 'ineq', 'fun': lambda x: e - d @ x}
bd = [(0, np.inf) for i in range(c1.shape[1])]
res = minimize(obj, np.random.randn(c1.shape[1]), constraints=constraints, bounds=bd)
print('最优解为：',np.round(res.x,4))
print('最优值为：',np.round(res.fun,4))
例 5.6 求下列非线性规划
minf(x)=x12+x22+x32+8,\min f(\boldsymbol{x}) = x_1^2 + x_2^2 + x_3^2 + 8,minf(x)=x12+x22+x32+8,
s. t. {x12−x2+x32≥0,x1+x22+x33≤20,−x1−x22+2=0,x2+2x32=3,x1,x2,x3≥0.\left\{\begin{array}{l}x_1^2 - x_2 + x_3^2 \geq 0, \\x_1 + x_2^2 + x_3^3 \leq 20, \\-x_1 - x_2^2 + 2 = 0, \\x_2 + 2x_3^2 = 3, \\x_1, x_2, x_3 \geq 0.\end{array}\right.⎩⎨⎧x12−x2+x32≥0,x1+x22+x33≤20,−x1−x22+2=0,x2+2x32=3,x1,x2,x3≥0.
解 求得当 x1=0.5522,x2=1.2033,x3=0.9478x_1 = 0.5522, x_2 = 1.2033, x_3 = 0.9478x1=0.5522,x2=1.2033,x3=0.9478 时，最小值 y=10.6511y = 10.6511y=10.6511。
from scipy.optimize import minimize
import numpy as np# 定义目标函数
obj = lambda x: sum(x ** 2) + 8# 不等式约束函数
def constraint1(x):x1, x2, x3 = xreturn [x1 ** 2 - x2 + x3 ** 2,20 - x1 - x2 ** 2 - x3 ** 3]# 等式约束函数
def constraint2(x):x1, x2, x3 = xreturn [-x1 - x2 ** 2 + 2,x2 + 2 * x3 ** 2 - 3]# 定义约束条件
constraints1 = {'type': 'ineq', 'fun': constraint1}
constraints2 = {'type': 'eq', 'fun': constraint2}
constraints = [constraints1, constraints2]
# 边界值
bd = [(0,np.inf) for i in range(3)]
res = minimize(obj, np.random.randn(3), constraints=constraints, bounds=bd)
print('最优解为：', np.round(res.x, 4))
print('最优值为：', np.round(res.fun, 4))
例 5.7 求解下列规划问题
minz=∣x1∣+2∣x2∣+3∣x3∣+4∣x4∣,s.t.{x1−x2−x3+x4=0,x1−x2+x3−3x4=1,x1−x2−2x3+3x4=−12.
\begin{aligned}
& \min z = |x_1| + 2|x_2| + 3|x_3| + 4|x_4|, \\
& \text{s.t.} \begin{cases}
x_1 - x_2 - x_3 + x_4 = 0, \\
x_1 - x_2 + x_3 - 3x_4 = 1, \\
x_1 - x_2 - 2x_3 + 3x_4 = -\frac{1}{2}.
\end{cases}
\end{aligned}
minz=∣x1∣+2∣x2∣+3∣x3∣+4∣x4∣,s.t.⎩⎨⎧x1−x2−x3+x4=0,x1−x2+x3−3x4=1,x1−x2−2x3+3x4=−21.
上述非线性规划问题是一个凸规划，可以使用 cvxpy 库求解，求得的最优解为 $ x_1 = 0.25 $, $ x_2 = x_3 = 0 $, $ x_4 = -0.25 $，目标函数的最优值为 1.25。
计算的 Python 程序如下：
import cvxpy as cp
import numpy as np# 定义变量
x = cp.Variable(4)# 定义目标向量
c = np.arange(1,5)
a = np.array([[1,-1,-1,1],[1,-1,1,-3],[1,-1,-2,3]])
b = np.array([0,1,-1/2])# 定义目标函数
obj = cp.Minimize(c @ cp.abs(x))
# 定义约束条件
constraints = [a @ x == b]
prob = cp.Problem(obj, constraints)
prob.solve(solver='GLPK_MI')
print('最优解为：',x.value)
print('最优值为：',prob.value)
例 5.8 (供应与选址) 建筑工地的位置（用平面坐标 a,ba,ba,b 表示，距离单位：km）及水泥日用量 ccc （单位：t）由表 5.4 给出。拟建两个料场向 6 个工地运送水泥，两个料场日储量各为 20t，问料场建在何处，使总的吨公里数最小。
表 5.4 建筑工地的位置及水泥日用量表
| 参数 | 工地 1 | 工地 2 | 工地 3 | 工地 4 | 工地 5 | 工地 6 |
|---|---|---|---|---|---|---|
| a/kma / \mathrm{km}a/km | 1.25 | 8.75 | 0.5 | 3.75 | 3 | 7.25 |
| b/kmb / \mathrm{km}b/km | 1.25 | 0.75 | 4.75 | 5 | 6.5 | 7.75 |
| c/tc / \mathrm{t}c/t | 3 | 5 | 4 | 7 | 6 | 11 |
解 记第 iii 个工地的位置为 (ai,bi)(i=1,2,⋯ ,6)(a_i, b_i) (i=1,2,\cdots,6)(ai,bi)(i=1,2,⋯,6)，水泥日用量为 cic_ici；拟建料场位置为 (xj,yj)(j=1,2)(x_j, y_j) (j=1,2)(xj,yj)(j=1,2)，日储量为 eje_jej，从料场 jjj 向工地 iii 的运送量为 zijz_{ij}zij。
建立如下的非线性规划模型：
min∑i=16∑j=12zij(xj−ai)2+(yj−bi)2,s.t.{∑j=12zij=ci,i=1,2,⋯ ,6,∑i=16zij≤ej,j=1,2,zij≥0,i=1,2,⋯ ,6;j=1,2.
\begin{aligned}
& \min \sum_{i=1}^{6} \sum_{j=1}^{2} z_{ij} \sqrt{(x_j - a_i)^2 + (y_j - b_i)^2}, \\
& \text{s.t.} \begin{cases}
\sum_{j=1}^{2} z_{ij} = c_i, & i = 1, 2, \cdots, 6, \\
\sum_{i=1}^{6} z_{ij} \leq e_j, & j = 1, 2, \\
z_{ij} \geq 0, & i = 1, 2, \cdots, 6; j = 1, 2.
\end{cases}
\end{aligned}
mini=1∑6j=1∑2zij(xj−ai)2+(yj−bi)2,s.t.⎩⎨⎧∑j=12zij=ci,∑i=16zij≤ej,zij≥0,i=1,2,⋯,6,j=1,2,i=1,2,⋯,6;j=1,2.
利用 Python 软件，求得拟建料场的坐标为 $ (3.2654, 5.1919) ，，， (7.25, 7.75) $。由两个料场向6个工地运料方案如表5.5所列,总的吨公里数为71.9352。
表 5.5 两个料场向6个工地运料方案
| 料 场 | 工地 1 | 工地 2 | 工地 3 | 工地 4 | 工地 5 | 工地 6 |
|---|---|---|---|---|---|---|
| 料场 1 | 0 | 5 | 0 | 0 | 0 | 11 |
| 料场 2 | 3 | 0 | 4 | 7 | 6 | 0 |
from scipy.optimize import minimize
import numpy as np# 读数据
data = np.loadtxt('建筑工地的位置及水泥日用量表.txt')
a = data[0];b = data[1];c = data[2]
e = np.array([20,20])# 定义目标函数
def obj(xyz):x = xyz[:2];y = xyz[2:4]z = xyz[4:].reshape(6, 2)obj = 0for i in range(6):for j in range(2):obj += z[i, j] * np.sqrt((x[j] - a[i]) ** 2 + (y[j] - b[i]) ** 2)return obj# 定义约束条件
construction = [{'type': 'eq', 'fun': lambda z: z[4:].reshape(6, 2).sum(axis=1) - c},{'type': 'ineq', 'fun': lambda z: e - z[4:].reshape(6, 2).sum(axis=0)}]# 边界
bd = [(0, np.inf) for i in range(16)]
res = minimize(obj, np.random.randn(16), constraints=construction, bounds=bd)
print('最优解为：', np.round(res.x, 4))
print('最优值为：', np.round(res.fun, 4))
