在Python中用nlopt添加矩阵系数多约束时遇索引越界问题求助

你遇到的IndexError核心原因是约束函数与nlopt的约束添加方法不匹配：

• 你编写的constraints函数是用于生成多个约束条件的（填充result数组），但你使用的add_inequality_constraint仅支持单个标量约束。这种 mismatch 导致nlopt传递了一个空的result数组，当你尝试向其写入值时就会触发索引越界错误。
• 你需要使用add_inequality_mconstraint（"m"代表multiple）来定义多变量约束，同时必须提前告知nlopt约束条件的总数，让它初始化正确大小的result数组。

import numpy as np
import sympy as sp
import nlopt

Qgreen = np.array([[1, 2, 3], 
                  [4, 5, 6], 
                  [7, 8, 9]])

len1, len2 = Qgreen.shape

# Calculate strength along the rows
strength1 = [np.sum(Qgreen[i, :]) for i in range(len1)]

# Calculate strength along the columns
strength2 = [np.sum(Qgreen[:, i]) for i in range(len2)]

# Combine the strengths into one list
str_combined = strength1 + strength2

# Define symbolic variables
vars_sym = sp.symbols(f'x1:{len(str_combined)+1}')  # x1, x2, ..., x6（重命名避免和内置vars冲突）

# Create the function using symbolic expressions
f_sum = sum(s * sp.log(v) for s, v in zip(str_combined, vars_sym))

# Flatten the product terms and their logarithms
product_terms = []
for i in range(len1):
    for j in range(len1, len1 + len2):
        product_terms.append(vars_sym[i] * vars_sym[j])

# Adding the log(1 - product_terms) part of the expression
f_product = sum(sp.log(1 - term) for term in product_terms)

# Final function combining both parts
f = f_sum + f_product

# Create a lambda function for fast evaluation
f_lambdified = sp.lambdify(vars_sym, f, "numpy")

# Define the function for nlopt (minimization)
def nlopt_function(x, grad):
    if grad.size > 0:
        # 预生成梯度的lambda函数提升性能（替代原循环evalf）
        grad_funcs = [sp.lambdify(vars_sym, sp.diff(f, var), "numpy") for var in vars_sym]
        grad[:] = np.array([func(*x) for func in grad_funcs])
    return f_lambdified(*x)

def constraints(x, result):
    count = 0
    num_vars = len(x)
    # 生成所有i≠j的约束（可选优化：只保留i<j避免重复约束）
    for i in range(num_vars):
        for j in range(num_vars):
            if i != j:
                result[count] = x[i] * x[j] - 1  # 约束：x[i]*x[j] < 1 → x[i]*x[j]-1 ≤ 0
                count += 1

# 计算变量总数和约束总数
num_vars = len(vars_sym)
num_constraints = num_vars * (num_vars - 1)  # 所有i≠j的组合

opt = nlopt.opt(nlopt.LN_COBYLA, num_vars)  # COBYLA支持不等式约束
opt.set_max_objective(nlopt_function)
# 使用多约束方法，指定约束数量和容忍度
opt.add_inequality_mconstraint(constraints, np.full(num_constraints, 1e-8))

initial_guess = np.random.uniform(0, 1, num_vars)
optimized_vars = opt.optimize(initial_guess)
max_val = opt.last_optimum_value()

print("Optimized variables:", optimized_vars)
print("Maximum value:", max_val)

1. 切换到多约束方法：

   • 用add_inequality_mconstraint替代add_inequality_constraint，明确告知nlopt约束条件的总数。
   • 传递与约束数量匹配的容忍度数组（或单个值自动广播）。

2. 优化梯度计算：

   • 原代码中循环调用evalf的方式性能较差，改为预生成梯度的lambda函数，提升优化速度。

3. 可选约束优化：

   • 原代码生成了重复约束（如(i,j)和(j,i)），因为x[i]*x[j]与x[j]*x[i]等价。若想减少约束数量，可将内层循环改为range(i+1, num_vars)，同时将num_constraints改为num_vars*(num_vars-1)//2。

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