Scipy带约束与边界的minimize：一般均衡模型优化问题

一、修正均衡条件与优化框架的耦合逻辑 #

不要在equilibrium函数内部做L_m的归一化操作，而是把所有约束显式地加入优化问题：

• 等式约束1：np.sum(L_m) == 1
• 等式约束2：np.multiply(G, L_m) - np.multiply(B, L_m) - C == 0（逐元素等于0）
• 不等式约束：L_m > 0（可用L_m >= 1e-8避免数值计算中的零值问题）

核心逻辑：equilibrium函数应直接输出由均衡方程求解得到的未归一化L_m，将均衡方程本身作为约束交给优化器处理，而非在函数内强行调整L_m破坏原均衡关系。

二、优化trust-constr的约束定义与参数设置 #

1. 正确定义非线性约束 #

使用scipy.optimize.NonlinearConstraint定义均衡相关约束（假设G、B依赖于x1，属于非线性约束）：

from scipy.optimize import NonlinearConstraint

# 均衡方程约束：确保G*L_m - B*L_m - C逐元素为0
def equilibrium_constraint(x1):
    L_m = equilibrium(x1)
    return np.multiply(G, L_m) - np.multiply(B, L_m) - C

eq_constraint1 = NonlinearConstraint(equilibrium_constraint, lb=0.0, ub=0.0, tol=1e-6)

# 劳动力份额总和约束
def sum_constraint(x1):
    L_m = equilibrium(x1)
    return np.sum(L_m) - 1

eq_constraint2 = NonlinearConstraint(sum_constraint, lb=0.0, ub=0.0, tol=1e-8)

# 劳动力份额非负约束
def pos_constraint(x1):
    L_m = equilibrium(x1)
    return L_m

ineq_constraint = NonlinearConstraint(pos_constraint, lb=1e-8, ub=np.inf, tol=1e-8)

2. 收紧收敛容差参数 #

默认容差可能过于宽松，导致算法提前终止。调用minimize时设置更严格的参数：

from scipy.optimize import minimize

result = minimize(obj1, x0, method='trust-constr',
                  constraints=[eq_constraint1, eq_constraint2, ineq_constraint],
                  options={
                      'gtol': 1e-8,  # 梯度收敛容差
                      'xtol': 1e-8,  # 变量变化容差
                      'barrier_tol': 1e-8,  # 障碍函数容差
                      'maxiter': 10000  # 增加最大迭代次数避免提前停止
                  })

3. 提供精确梯度/雅可比矩阵 #

trust-constr算法在有精确梯度时收敛性会大幅提升。若能推导解析梯度，务必提供；无法推导时，需确保数值梯度的步长足够小：

# 示例：为目标函数和约束定义解析梯度/雅可比
def obj1_grad(x1):
    # 返回与x1同维度的目标函数梯度向量
    pass

def equilibrium_constraint_jac(x1):
    # 返回均衡约束的雅可比矩阵，形状为(约束维度, x1维度)
    pass

# 在约束中指定jac参数
eq_constraint1 = NonlinearConstraint(equilibrium_constraint, lb=0.0, ub=0.0,
                                     jac=equilibrium_constraint_jac, tol=1e-6)

# 调用minimize时传入目标函数梯度
result = minimize(obj1, x0, method='trust-constr', jac=obj1_grad,
                  constraints=[eq_constraint1, eq_constraint2, ineq_constraint],
                  options={'gtol': 1e-8, 'xtol': 1e-8, 'maxiter': 10000})

三、验证优化结果的约束满足情况 #

优化完成后，手动验证约束是否达标：

L_m_opt = equilibrium(result.x)
print("均衡方程残差：", np.linalg.norm(np.multiply(G, L_m_opt) - np.multiply(B, L_m_opt) - C))
print("劳动力份额总和：", np.sum(L_m_opt))
print("劳动力份额均为正：", np.all(L_m_opt > 1e-8))

若残差仍较大，需检查equilibrium函数的实现逻辑——确认L_m是否是从均衡方程中正确解出的，而非近似值。

四、备选方案：带惩罚项的目标函数重构 #

若trust-constr仍不收敛，可将约束嵌入目标函数使用惩罚法，再用L-BFGS-B求解：

def obj_with_penalty(x1, lambda1=1e6, lambda2=1e6):
    L_m = equilibrium(x1)
    # 原目标函数值
    obj_val = obj1(x1)
    # 均衡方程的平方惩罚项
    penalty1 = lambda1 * np.linalg.norm(np.multiply(G, L_m) - np.multiply(B, L_m) - C)**2
    # 份额和为1的平方惩罚项
    penalty2 = lambda2 * (np.sum(L_m) - 1)**2
    # 非负约束的对数障碍惩罚
    penalty3 = -lambda1 * np.sum(np.log(L_m + 1e-8))
    return obj_val + penalty1 + penalty2 + penalty3

注意惩罚系数需逐步调整，从较小值开始递增，避免数值不稳定。

内容的提问来源于stack exchange，提问作者Alien