一般的非線性規劃求解(非凸函數)

例 5.5 (續例 5.3) 求解如下二次規劃模型
$$\begin{array}{rl}& \min ?x_1^2?0.3x_1{x}_2?2x_2^2+98x_1+277x_2\\ & \text{s.t.}?\begin{array}{l}{x}_1+x_2\le 100,\\ x_1?2x_2\le 0,\\ x_1,x_2\ge 0.\end{array}\end{array}$$
由于上述二次規劃的目標函數不是凸函數，不能使用 cvxpy 庫進行求解。

使用 Python 軟件求得的最優解為
$$x_1=0,x_2=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 求下列非線性規劃

$\min f(\mathbf{x})=x_1^2+x_2^2+x_3^2+8,$

s. t. $?\begin{array}{l}{x}_1^2?x_2+x_3^2\ge 0,\\ x_1+x_2^2+x_3^3\le 20,\\ ?x_1?x_2^2+2=0,\\ x_2+2x_3^2=3,\\ x_1,x_2,x_3\ge 0.\end{array}$

解 求得當 $x_1=0.5522,x_2=1.2033,x_3=0.9478$ 時，最小值 $y=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 求解下列規劃問題
$$\begin{array}{rl}& \min z=|x_1|+2|x_2|+3|x_3|+4|x_4|,\\ & \text{s.t.}?\begin{array}{l}{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{array}\end{array}$$

上述非線性規劃問題是一個凸規劃，可以使用 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,b$ 表示，距離單位：km）及水泥日用量 $c$ （單位：t）由表 5.4 給出。擬建兩個料場向 6 個工地運送水泥，兩個料場日儲量各為 20t，問料場建在何處，使總的噸公里數最小。

表 5.4 建筑工地的位置及水泥日用量表

參數	工地 1	工地 2	工地 3	工地 4	工地 5	工地 6
$a/\mathrm{km}$	1.25	8.75	0.5	3.75	3	7.25
$b/\mathrm{km}$	1.25	0.75	4.75	5	6.5	7.75
$c/\mathrm{t}$	3	5	4	7	6	11

解 記第 $i$ 個工地的位置為 $(a_i,b_i)(i=1,2,?,6)$，水泥日用量為 $c_i$；擬建料場位置為 $(x_j,y_j)(j=1,2)$，日儲量為 $e_j$，從料場 $j$ 向工地 $i$ 的運送量為 $z_{ij}$。

建立如下的非線性規劃模型：
$$\begin{array}{rl}& \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{array}{ll}{\sum }_{j=1}^2{z}_{ij}=c_i,& i=1,2,?,6,\\ {\sum }_{i=1}^6{z}_{ij}\le e_j,& j=1,2,\\ z_{ij}\ge 0,& i=1,2,?,6;j=1,2.\end{array}\end{array}$$
利用 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))