带约束整数非线性优化报错：zero-size array to reduction operation maximum

我正在求解一个带两个约束条件的整数非线性优化问题，运行代码时触发如下错误：

ValueError: zero-size array to reduction operation maximum which has no identity

以下是完整代码结构：

导入依赖与初始化数据 #

from ypstruct import structure
import numpy as np
import xlrd
from scipy.optimize import minimize

# 输入数据
input_data= np.array([[0.009, 0.001, 0.001, 0.009, 0.001, 0.009, 0.001, 0.009, 0.001, 0.001, 0.0002, 0.0002, 0.001, 0.001, 0.001, 0.009, 0.005, 0.001],
                     [0.05, 0.05, 0.12, 0.12, 0.15, 0.15, 0.3, 0.3, 0.02, 0.02, 0.18, 0.18, 0.2, 0.2, 0.01, 0.01, 0.01, 0.01],
                     [60, 40, 60, 40, 60, 40, 60, 40, 60, 40, 60, 40, 60, 40, 60, 40, 60, 40],   
                     [30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30],
                     [2000,2000,1992,1992,1965,1965,1976,1976,1979,1979,1982,1982,1984,1984,1994,1994,1991,2006],
                     [5,5,3,3,8,8,8,8,10,10,12,12,11,11,1,1,1,1]])
ngroup = 9 ; limit_per_year=9 ; starting_year = 2000

C_capital = np.zeros(ngroup*2); C_sudden_replace = np.zeros(ngroup*2)
C_Gain_Power = np.zeros(ngroup*2); penalty=np.zeros(ngroup*2)
finalcost =np.zeros(ngroup*2)
X= np.array([32,31,36,38,49,44,54,54,44,45,34,34,24,24,38,85,45,54])

目标函数 #

def objective (X):
   discount_rate= 0.04; TOTAL_COST=0
   penalti= 5
   for i in range(ngroup*2):
       C_capital[i]= (X[i]/((1+X[i])**discount_rate))
       C_sudden_replace[i] =(np.exp(-input_data[1,i]*X[i]))*input_data[2,i]
       C_Gain_Power[i] = (50 - X[i])* input_data[1,i]
       penalty[i] = 0
       if i % 2 ==1:
           penalty[i]= penalti*(abs((X[i]+input_data[4,i]-(X[i-1]+input_data[4,i]))))

       finalcost[i] = C_capital[i] + C_sudden_replace[i] + C_Gain_Power[i] + penalty[i]
       TOTAL_COST=TOTAL_COST+ finalcost[i]

   return TOTAL_COST

COSTOBJ = objective(X)
print(COSTOBJ)

第一个约束条件 #

def constraint1(X):
    matrix1= np.zeros([ngroup*2,12])
    for i in range(ngroup*2):
        d=0
        for j in range(12):
            if j<= input_data[5,i]:
               matrix1[i,j]=input_data[4,i]+X[i]+d
               d=d+1
    #print(matrix1)
    a1=int(np.min(matrix1[np.nonzero(matrix1)]))
    a2=int(np.max(matrix1[np.nonzero(matrix1)]))
    matrix2=np.zeros((1 ,int(a2-a1)+1))
    #print(matrix2)
    b=0
    for i in range(a1, a2+1):
        matrix2[0,b]=np.count_nonzero(matrix1 == i)
        b=b+1
    sortir1= int(np.max(matrix2[np.nonzero(matrix2)]))
    return limit_per_year-sortir1

CONS1= constraint1(X)
print(CONS1)

第二个约束条件 #

def constraint2(X):
    matrix1= np.zeros([ngroup*2,12])
    for i in range(ngroup*2):
        d=0
        for j in range(12):
            if j<= input_data[5,i]:
               matrix1[i,j]=input_data[4,i]+X[i]+d
               d=d+1
    #print(matrix1)
    sortir= int(np.min(matrix1[:, 0]))
    return sortir-starting_year

CONS2= constraint2(X)
print(CONS2)

优化执行代码（触发错误） #

cons1 = {'type': 'ineq', 'fun': constraint1}
cons2 = {'type': 'ineq', 'fun': constraint2}

cons=[cons1,cons2]
sol= minimize(objective,X ,method='SLSQP', constraints=cons)

错误原因分析 #

这个错误本质是调用np.max()或np.min()时传入了空数组。具体到你的代码：

• 在constraint1中，当优化器尝试某些X值时，可能导致matrix1中没有非零元素，此时matrix1[np.nonzero(matrix1)]返回空数组，np.min()和np.max()就会抛出这个错误。
• 另外，constraint2中直接取matrix1[:,0]的最小值，若该列全为0，虽然不会触发错误，但可能不符合你的约束逻辑。

解决方法 #

1. 给空数组的极值计算加兜底逻辑

修改constraint1，处理空数组的情况：

def constraint1(X):
    matrix1= np.zeros([ngroup*2,12])
    for i in range(ngroup*2):
        d=0
        for j in range(12):
            if j<= input_data[5,i]:
               matrix1[i,j]=input_data[4,i]+X[i]+d
               d=d+1
    # 提取非零元素，判断是否为空
    non_zero_matrix1 = matrix1[np.nonzero(matrix1)]
    if len(non_zero_matrix1) == 0:
        # 根据业务逻辑返回一个满足约束的值，比如直接返回limit_per_year
        return limit_per_year
    
    a1=int(np.min(non_zero_matrix1))
    a2=int(np.max(non_zero_matrix1))
    matrix2=np.zeros((1 ,int(a2-a1)+1))
    b=0
    for i in range(a1, a2+1):
        matrix2[0,b]=np.count_nonzero(matrix1 == i)
        b=b+1
    
    # 同样处理matrix2的空数组情况
    non_zero_matrix2 = matrix2[np.nonzero(matrix2)]
    sortir1= int(np.max(non_zero_matrix2)) if len(non_zero_matrix2) >0 else 0
    return limit_per_year-sortir1

2. 给X添加边界约束，避免异常值

优化器可能会尝试不合理的X值（比如负数），导致input_data[4,i]+X[i]+d计算结果为0或负数，进而让matrix1出现全零行/列。给X设置合理的取值范围：

# 假设每个X元素的取值范围是1到100，可根据实际业务调整
bounds = [(1, 100) for _ in range(len(X))]
sol= minimize(objective,X ,method='SLSQP', constraints=cons, bounds=bounds)

3. 注意整数优化的局限性

scipy.optimize.minimize的SLSQP方法是针对连续变量的优化，如果你需要严格的整数解，直接运行得到的结果是连续值，需要后续取整处理，或者改用专门的整数优化工具（比如pyscipopt、Google OR-Tools）。

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