SciPy 的 linprog 函数是求解线性规划问题的高效工具,它能处理线性目标函数以及由矩阵和向量表达的一组线性等式或不等式约束。然而在实际建模中,很多业务限制并不是线性的,例如风险控制中的方差上限、物流网络里的距离半径限制、设备运行中的功率平方和约束等。此时如果仍然希望沿用线性目标函数,只是把约束条件扩展到非线性形式,就需要离开 linprog,转而使用 SciPy 更通用的优化接口 scipy.optimize.minimize。该接口支持多种约束类型,其中 NonlinearConstraint 可以方便地封装自定义非线性函数,并设置约束上下界。

一、线性规划求解器与非线性的边界
标准线性规划问题通常写作最小化 c^T x,并满足 A_ub x ≤ b_ub、A_eq x = b_eq 以及变量边界 lb ≤ x ≤ ub。SciPy 中的 linprog 默认使用 HiGHS 求解器,能够快速处理大规模稀疏线性模型。比如下面这个简单例子,目标是最小化 -x0 + 4x1,同时满足 -3x0 + x1 ≤ 6 和 x0 + 2x1 ≤ 4,并且 x1 的下界为 -3。
from scipy.optimize import linprog c = [-1, 4] A_ub = [[-3, 1], [1, 2]] b_ub = [6, 4] bounds = [(None, None), (-3, None)] res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs') print(res.x) print(res.fun)
这段代码可以顺利求解,因为所有表达式都是决策变量的一次函数。一旦约束中出现 x0 的平方、x0 乘 x1 的交互项,或者 sin、exp 等非线性函数,线性规划模型就不再适用。比如把约束改成 x0^2 + x1^2 ≤ 1,试图用 linprog 去求解就会报错,因为它无法理解这种约束形式。这个边界要求开发者在建模阶段就识别出哪些限制属于非线性约束,并选择对应的求解策略。
一种常见错误是试图把非线性约束强行线性化,例如把圆形区域 x0^2 + x1^2 ≤ 1 改成四个线性不等式 -1 ≤ x0 ≤ 1、-1 ≤ x1 ≤ 1,虽然这样能求解,但可行域被明显放大,得到的结果可能完全不可用。正确做法是保留非线性约束,并使用支持非线性约束的求解器。
二、使用 minimize 与 NonlinearConstraint 接入自定义约束
SciPy 的 scipy.optimize.minimize 提供了 SLSQP、trust-constr、COBYLA 等支持约束优化的方法。其中 SLSQP 是最常用的序列二次规划方法,适合中小规模问题;trust-constr 则适合带稀疏雅可比矩阵的较大规模问题。为了把非线性约束描述清楚,可以使用 NonlinearConstraint 对象。它的第一个参数是约束函数,该函数接收决策变量数组,返回一个标量或数组;第二个和第三个参数分别是约束的下界和上界。
下面给出一个例子:目标函数仍然保持线性,为 -2x0 - 3x1,等价于最大化 2x0 + 3x1;非线性约束为 x0^2 + x1^2 ≤ 1,表示决策点必须落在单位圆内;同时还有一个线性约束 x0 + x1 ≥ 0.5。使用 SLSQP 求解时代码如下。
import numpy as np
from scipy.optimize import minimize, NonlinearConstraint, LinearConstraint
def objective(x):
return -2 * x[0] - 3 * x[1]
def circle_constraint(x):
return x[0] ** 2 + x[1] ** 2
nlc = NonlinearConstraint(circle_constraint, -np.inf, 1)
A = [[1, 1]]
lb_linear = [0.5]
linear_constraint = LinearConstraint(A, lb_linear, np.inf)
x0 = [0.6, 0.6]
res = minimize(objective, x0, method='SLSQP',
constraints=[nlc, linear_constraint],
bounds=[(-2, 2), (-2, 2)])
print(res.x)
print(res.fun)
运行这段代码后,求解器会尝试在单位圆内且满足 x0 + x1 ≥ 0.5 的区域内找到使得 -2x0 - 3x1 最小的点。由于目标函数是线性的,最优解通常会落在约束边界上。这里 NonlinearConstraint 的上界设为 1,下界设为负无穷,表示只要求圆约束函数值不超过 1。如果需要约束函数值恰好等于某个常数,可以把下界和上界设为同一个数。
相比 linprog,minimize 的优势是约束定义方式更灵活,可以在同一个问题里混合线性约束、非线性约束以及变量边界。不过它也要求用户提供一个初始点 x0,而初始点的选择会直接影响求解效率和最终结果。
三、约束函数与雅可比矩阵的编写方式
对于非线性约束,求解器默认使用有限差分方法估计梯度,这在变量数量较少时足够稳定。但当变量维度较大,或者约束函数计算代价很高时,有限差分会显著增加运行时间。此时可以提供解析雅可比矩阵,帮助求解器更快收敛。雅可比矩阵的每一行对应一个约束,列对应决策变量,元素为约束函数对该变量的偏导数。
以圆形约束 x0^2 + x1^2 为例,它只返回一个标量,因此雅可比矩阵的形状是 1×2,元素分别为 2x0 和 2x1。可以像下面这样定义并传入 NonlinearConstraint。
def circle_constraint(x):
return np.array([x[0] ** 2 + x[1] ** 2])
def circle_jac(x):
return np.array([[2 * x[0], 2 * x[1]]])
nlc = NonlinearConstraint(circle_constraint, -np.inf, 1, jac=circle_jac)
如果同时有多个非线性约束,约束函数可以返回一个数组。例如增加约束 exp(x0) - x1 ≤ 2,那么约束函数返回包含两项的数组,雅可比矩阵的行数也要对应增加。雅可比矩阵的第一行是圆形约束的偏导数,第二行是 exp(x0) - x1 的偏导数,即 [exp(x0), -1]。
使用 trust-constr 方法时,还可以提供 Hessian 矩阵来进一步加速收敛,不过对于大多数中小规模问题,提供雅可比矩阵已经足够。需要注意的是,解析雅可比矩阵里的每个元素都要与约束函数的输出顺序严格对应,否则求解器可能出现收敛困难或结果错误。
四、初始点、可行性与局部最优风险
非线性约束会使可行域变成非凸集合,比如两个圆形可行域的交集可能不再是一个凸集。SLSQP 这类局部优化算法在这样的非凸可行域上容易陷入局部最优,而无法保证找到全局最优解。一个常见的规避策略是多起点搜索:在可行域内生成多个初始点,分别调用 minimize,然后选择目标函数值最优的那个结果。虽然这样不能理论上保证全局最优,但在工程实践中往往足够可靠。
初始点的选择也很关键。SLSQP 在迭代初期会尝试修复约束违反,但如果初始点离可行域太远,求解器可能无法收敛。最好根据业务经验给出一个大致可行的初始点。例如在资源分配问题中,可以先把资源平均分配,确保满足总量约束,再检查非线性约束是否满足。如果不满足,可以适当缩小分配比例,直到进入可行域。
如果问题规模不大,但非线性约束较多且可行域复杂,还可以考虑使用 differential_evolution 或 shgo 等全局优化方法。它们不依赖梯度,能够探索更大的搜索空间,但计算成本通常远高于 SLSQP。因此建议先用 SLSQP 快速求解,结合多起点策略观察结果稳定性,再决定是否需要引入全局优化器。
五、圆形约束下的资源分配实例
下面给出一个更贴近实际场景的例子。假设有三种资产配置,决策变量 x0、x1、x2 分别表示投入比例,三者之和必须等于 1。线性目标是最小化负收益,等价于最大化收益 -0.10x0 - 0.15x1 - 0.08x2。此外还有一个非线性风险约束,要求组合方差不超过 0.02,方差函数为 0.04x0^2 + 0.09x1^2 + 0.01x2^2 + 0.02x0x1。这个约束是二次的,无法放进 linprog,但可以通过 SLSQP 处理。
import numpy as np
from scipy.optimize import minimize, NonlinearConstraint, LinearConstraint
def objective(x):
return -(0.10 * x[0] + 0.15 * x[1] + 0.08 * x[2])
def variance_constraint(x):
return (0.04 * x[0] ** 2 + 0.09 * x[1] ** 2
+ 0.01 * x[2] ** 2 + 0.02 * x[0] * x[1])
def variance_jac(x):
return np.array([[0.08 * x[0] + 0.02 * x[1],
0.18 * x[1] + 0.02 * x[0],
0.02 * x[2]]])
nlc = NonlinearConstraint(variance_constraint, -np.inf, 0.02,
jac=variance_jac)
A = [[1, 1, 1]]
linear_constraint = LinearConstraint(A, [1], [1])
x0 = [0.33, 0.33, 0.34]
res = minimize(objective, x0, method='SLSQP',
constraints=[nlc, linear_constraint],
bounds=[(0, 1), (0, 1), (0, 1)])
print(res.x)
print('最大收益:', -res.fun)
print('组合方差:', variance_constraint(res.x))
这段代码把线性等式约束 x0 + x1 + x2 = 1 通过 LinearConstraint 实现,把风险上限通过 NonlinearConstraint 实现。求解结果会显示三种资产的最优配置比例、最大收益以及对应的组合方差。由于目标函数是线性的,而风险约束是凸的二次函数,该问题仍然是一个凸优化问题,SLSQP 能够稳定地找到全局最优解。
从求解结果还可以观察到一个重要现象:当风险上限收紧到一定程度时,最优解会落在风险约束的边界上;当风险上限很宽松时,最优解会落在变量边界或线性等式约束的边界上。这说明非线性约束是否起作用,取决于它是否成为限制目标改进的活跃约束。理解这一点有助于在模型调试阶段判断约束设置是否合理。