Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

2.2 线性规划

线性规划是指目标函数和所有约束条件均为变量的线性函数的优化问题。

案例：Chebyshev中心问题

本节通过寻找Chebyshev中心的例子来展示线性规划的应用。

Chebyshev中心问题旨在寻找一个半径最大，并且完全包含在给定多面体内部的球体。在二维平面中，它直观上等同于在一个多边形内部寻找“最大内切圆”。 这个问题可以极为精妙地通过将它的“圆心”和“半径”作为优化变量，转化为典型的线性规划问题进行求解。

cvxpy 版本

# 导入cvxpy 和 numpy
import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
x_grid = np.linspace(-1, 2, 10)
fig, ax = plt.subplots(figsize=(5,4))
ax.plot(x_grid, 5 - x_grid, label=r"$y=5-x$")
ax.plot(x_grid, 2 - x_grid, label=r"$y=2-x$")
ax.plot(x_grid, 2 + 1.5 *x_grid, label=r"$y=2+\frac{3}{2}x$")
ax.plot(x_grid, 1 + x_grid, label=r"$y=1+x$")
ax.legend(loc='upper left',fontsize=12)
ax.set_aspect('equal')
ax.grid(True)
ax.set_xlim([-1,2])
ax.set_ylim([1,4])
plt.tight_layout()
<Figure size 500x400 with 1 Axes>
a1 = np.array([1,1])
b1 = 5
a2 = np.array([-1,-1])
b2 = -2
a3 = np.array([-1.5,1])
b3 = 2
a4 = np.array([1,-1])
b4 = -1
# 创建优化变量
x = cp.Variable(2)
r = cp.Variable(nonneg=True)

# 创建限制条件
constraints = [
    a1 @ x + r * np.linalg.norm(a1) <= b1,
    a2 @ x + r * np.linalg.norm(a2) <= b2,
    a3 @ x + r * np.linalg.norm(a3) <= b3,
    a4 @ x + r * np.linalg.norm(a4) <= b4
]

# 创建目标函数
obj = cp.Maximize(r)

# 创建优化问题
prob = cp.Problem(obj, constraints)
prob.solve()

# 输出结果
print("求解状态:", prob.status)
print("最优值 (半径 r):", f"{prob.value:.4f}" if prob.value is not None else None)
print("最优变量 圆心 x:", np.round(x.value, 4) if x.value is not None else None)
print("最优变量 半径 r:", np.round(r.value, 4) if r.value is not None else None)
求解状态: optimal
最优值 (半径 r): 0.5097
最优变量 圆心 x: [1.2792 3.    ]
最优变量 半径 r: 0.5097

fig, ax = plt.subplots(figsize=(5,4))
ax.plot(x_grid, 5 - x_grid, label=r"$y=5-x$")
ax.plot(x_grid, 2 - x_grid, label=r"$y=2-x$")
ax.plot(x_grid, 2 + 1.5 *x_grid, label=r"$y=2+\frac{3}{2}x$")
ax.plot(x_grid, 1 + x_grid, label=r"$y=1+x$")

circle = Circle(x.value, r.value, edgecolor='blue', facecolor='none', linewidth=2)

ax.add_patch(circle)
ax.legend(loc='upper left',fontsize=12)
ax.set_aspect('equal')
ax.grid(True)
ax.set_xlim([-1,2])
ax.set_ylim([1,4])
plt.tight_layout()
<Figure size 500x400 with 1 Axes>

scipy.optimize 版本

使用 scipy.optimize.linprog 求解 Chebyshev 中心问题。

与 CVXPY 的声明式建模不同，如果我们在底层使用 SciPy 的线性规划模块 linprog，我们需要将数学问题转化为标准形式：

min⁡cTx\min c^T \mathbf{x}
(1)

s.t. Aubx≤bub\text{s.t. } A_{ub} \mathbf{x} \le b_{ub}
(2)

在这个案例中，优化变量由坐标 (x1,x2)(x_1, x_2) 和半径 rr 组成，即 vars=[x1,x2,r]T\mathbf{vars} = [x_1, x_2, r]^T。

  • 目标函数 是最大化 rr，即最小化 −r-r。所以目标向量 c=[0,0,−1]Tc = [0, 0, -1]^T。

  • 约束条件：我们需要把每一个限制 aiTx+r∥ai∥≤bia_i^T x + r \|a_i\| \le b_i 拼成一个矩阵形式 Aubvars≤bubA_{ub} \mathbf{vars} \le b_{ub}。

from scipy.optimize import linprog

# 1. 定义目标函数向量 c，使得 c^T vars = -r (我们要最小化 -r)
# 变量顺序为: [x1, x2, r]
c = np.array([0, 0, -1])

# 2. 定义不等式约束 A_ub * vars <= b_ub
# 提示:约束形式为 a_i^T x + r * ||a_i|| <= b_i
# 对于每个约束，系数向量为 [a_i(1), a_i(2), ||a_i||]
A_ub = np.array([
    [a1[0], a1[1], np.linalg.norm(a1)],
    [a2[0], a2[1], np.linalg.norm(a2)],
    [a3[0], a3[1], np.linalg.norm(a3)],
    [a4[0], a4[1], np.linalg.norm(a4)]
])

b_ub = np.array([b1, b2, b3, b4])

# 3. 设置变量边界 (Bounds)
# 坐标 x1, x2 无限制 (负无穷到正无穷)
# 半径 r 必须大于等于 0
bounds = [(None, None), (None, None), (0, None)]

# 4. 调用 linprog 求解
res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds)

# 5. 提取并打印结果
print("求解状态 (成功):", res.success)
print("详细信息:", res.message)

if res.success:
    # 注意：我们的目标是最大化 r，而 linprog 计算的是最小化 -r
    # 因此真实的最优半径就是 -res.fun (或者直接从 res.x 提取 r)
    print("最优值 (半径 r):", f"{-res.fun:.4f}")
    print("最优变量 圆心 x:", np.round(res.x[0:2], 4))
    print("最优变量 半径 r:", np.round(res.x[2], 4))
求解状态 (成功): True
详细信息: Optimization terminated successfully. (HiGHS Status 7: Optimal)
最优值 (半径 r): 0.5097
最优变量 圆心 x: [1.2792 3.    ]
最优变量 半径 r: 0.5097

使用通用非线性求解器 scipy.optimize.minimize 取代线性规划专用版求解

为了拓展思路，虽然这是一个严格的线性规划问题，我们同样也可以“大材小用”，使用通用求解器 scipy.optimize.minimize来进行数值迭代求解。

💡 提示

对于纯线性的优化目标与约束体系，建议在生产环境优先选用基于单纯形法或内点法的线性专用求解器 (如 scipy.optimize.linprog，或其他专用 LP Solver)。

只有当你的环境只能暴露基础函数形式，或者约束中有非线性方程时，才退而求其次使用下面的通用数值求解器 minimize。

from scipy.optimize import minimize

# 1. 目标函数: 我们要最小化 -r (由于决策变量是 [x1, x2, r]，对应的索引为 [0, 1, 2])
def objective(vars):
    return -vars[2]

# 2. 约束函数 a_i * x + r * ||a_i|| <= b_i
# 对于 minimize，我们需要把不等式约束转换为 ">= 0" 的形式:
# 即: b_i - (a_i * x + r * ||a_i||) >= 0
def ineq_constraint(vars):
    x = vars[0:2]
    r = vars[2]
    
    # 也可以返回一个列表作为多个并行约束的输出
    return [
        b1 - (np.dot(a1, x) + r * np.linalg.norm(a1)),
        b2 - (np.dot(a2, x) + r * np.linalg.norm(a2)),
        b3 - (np.dot(a3, x) + r * np.linalg.norm(a3)),
        b4 - (np.dot(a4, x) + r * np.linalg.norm(a4))
    ]

# 包装成 solver 所需配置
# 也可以单独为半径 r 指定边界 bounds = [(None, None), (None, None), (0, None)]
constraints_minimize = {'type': 'ineq', 'fun': ineq_constraint}
bounds = [(None, None), (None, None), (0.0, None)]

# 3. 设定初始猜测点 (随便写一个位于空间内的点和正数r)
vars_0 = [1.0, 2.0, 0.5] 

# 4. 调用通用的 minimize 求解器
res_min = minimize(objective, vars_0, bounds=bounds, constraints=constraints_minimize)

# 输出结果
print("通用求解器求解状态 (Success):", res_min.success)
print("详细信息:", res_min.message)

if res_min.success:
    print("最优值 (半径 r):", f"{-res_min.fun:.4f}")
    print("最优变量 圆心 x:", np.round(res_min.x[0:2], 4))
    print("最优变量 半径 r:", np.round(res_min.x[2], 4))
通用求解器求解状态 (Success): True
详细信息: Optimization terminated successfully
最优值 (半径 r): 0.5097
最优变量 圆心 x: [1.2792 3.    ]
最优变量 半径 r: 0.5097
MATH3806 Code
2.1 CVXPY 和 SciPy 简介
MATH3806 Code
2.3 二次规划