如何修复优化代码中的RuntimeWarning与迭代超限AssertionError

运行优化问题代码时，遇到两个错误：

• RuntimeWarning: invalid value encountered in double_scalars
• AssertionError: Iteration limit reached
  曾怀疑是初始条件问题，但修改初始条件后错误依旧存在。

代码如下：

import numpy as np
import sympy as sp
from sympy import *
from scipy.special import roots_legendre, eval_legendre

from scipy.optimize import minimize, NonlinearConstraint
# create variables
N=8;n1=2;n2=3;n3=2;n4=3;z=1;m=2;a=0;b=1
x = [sp.symbols('x%d' % i) for i in range(2*(n1+n2+n3+n4)-4)]
s=zeros(2*n1-1,1);e=zeros(2*n2-1,1);h=zeros(2*n3-1,1);l=zeros(2*n4-1)
for i in range(2*n1-1):
    s[i]=x[i]
for i in range(2*n2-1):
    e[i]=x[i+2*n1-1]
for i in range(2*n3-1):
    h[i]=x[i+2*(n1+n2)-2]
for i in range(2*n4-1):
    l[i]=x[i+2*(n1+n2+n3)-3]   
#s = [sp.symbols('s%d' % i) for i in range(2*n1-1)]
#e = [sp.symbols('e%d' % i) for i in range(2*n2-1)]#Instead of "i" we choose "e"
#h = [sp.symbols('h%d' % i) for i in range(2*n3-1)]
#l = [sp.symbols('l%d' % i) for i in range(2*n4-1)]
from sympy.abc import t
B=zeros(N,N)
number = range(N)
for i in number:
   for j in number:
       if i<j :
           B[i,j]=0
       else :
           B[i,j]=sp.factorial(i+j)/(2**j*sp.factorial(j)*sp.factorial(i-j))
PS=Matrix(1,n1,s[0:n1]);PI=Matrix(1,n2,e[0:n2]);PH=Matrix(1,n3,h[0:n3]);PL=Matrix(1,n4,l[0:n4])
TS=[[1]];TI=[[1]];TH=[[1]];TL=[[1]]
for i in range(n1-1):
    TS.append([t**(s[n1+i]+i+1)])
for i in range(n2-1):    
    TI.append([t**(e[n2+i]+i+1)])
for i in range(n3-1):        
    TH.append([t**(h[n3+i]+i+1)])
for i in range(n4-1):            
    TL.append([t**(l[n4+i]+i+1)])
MTS=Matrix(TS);MTI=Matrix(TI);MTH=Matrix(TH);MTL=Matrix(TL)
BS=B[0:n1,0:n1];BI=B[0:n2,0:n2];BH=B[0:n3,0:n3];BL=B[0:n4,0:n4]
S=PS*BS*MTS;II=PI*BI*MTI;H=PH*BH*MTH;L=PL*BL*MTL

#CONVERT SYMPY MATRICES TO NUMPY ONE
S0=np.array(S);I0=np.array(II);H0=np.array(H);L0=np.array(L)
GS=zeros(n1,n1);GI=zeros(n2,n2);GH=zeros(n3,n3);GL=zeros(n4,n4)

for i in range(n1):
   if i==0:
      GS[i,i]=0
   else:
       GS[i,i]=simplify(gamma(s[n1+i-1]+1)/gamma(s[n1+i-1]+1-z))
for i in range(n2):
    if i==0:
      GI[i,i]=0
    else:
     GI[i,i]=simplify(gamma(e[n2+i-1]+1)/gamma(e[n2+i-1]+1-z))
for i in range(n3):
    if i==0:
      GH[i,i]=0
    else:
      GH[i,i]=simplify(gamma(h[n3+i-1]+1)/gamma(h[n3+i-1]+1-z))
for i in range(n4):
    if i==0:
      GL[i,i]=0
    else:
      GL[i,i]=simplify(gamma(l[n4+i-1]+1)/gamma(l[n4+i-1]+1-z))     
DS=simplify(t**(-z)*PS*BS*GS*MTS)
DI=simplify(t**(-z)*PI*BI*GI*MTI)
DH=simplify(t**(-z)*PH*BH*GH*MTH)
DL=simplify(t**(-z)*PL*BL*GL*MTL)
RS=DS[0,0]-(.0043217-.5944*S[0,0]*II[0,0]-(.025+3.5)*S[0,0])
RI=DI[0,0]-(.5944*S[0,0]*II[0,0]+.0056*H[0,0]*II[0,0]+.27*L[0,0]-(.025+3.5+.0025+.5)*II[0,0])
RH=DH[0,0]-(.535-.0056*H[0,0]*II[0,0]+.5*II[0,0]-(.025+3.5)*H[0,0])
RL=DL[0,0]-(.0025*II[0,0]-(.025+3.5+.27)*L[0,0])      
#f = lambdify(t, R)
import decimal
decimal.getcontext().prec=4
def R(k):
   return (RS**2+RI**2+RH**2+RL**2).subs(t,k)
r, w = roots_legendre(m)
F=zeros(1,m)
I1=0
for i in range (m):
   F[i]=r[i]*R((b-a)/2*r[i]+(b+a)/2)
   I1=I1+F[i]
OF1=((b-a)/2)*I1 #OF=Objective Function
OF=OF1.evalf(3)
t = next(s for s in S.free_symbols if s.name == 't')
x = [
    next(s for s in OF.free_symbols if s.name == f'x{i}')
    for i in range(2*(n1+n2+n3+n4)-4)
]
args = [t, *x]
OF = sp.lambdify(args, OF)
S = sp.lambdify(args, S.flat()[0])
I = sp.lambdify(args, II.flat()[0])
H = sp.lambdify(args, H.flat()[0])
L = sp.lambdify(args, L.flat()[0])

def objective(params: np.ndarray) -> float:
    return OF(*params)

def Sc(params: np.ndarray) -> float:
    return S(*params)
def Ic(params: np.ndarray) -> float:
    return I(*params)
def Hc(params: np.ndarray) -> float:
    return H(*params)
def Lc(params: np.ndarray) -> float:
    return L(*params)

res = minimize(
    fun=objective,
    x0=np.ones(1 + 2*(n1+n2+n3+n4)-4),
    constraints=(
        NonlinearConstraint(fun=Sc, lb=43994, ub=43994),
        NonlinearConstraint(fun=Ic, lb=1, ub=1),
        NonlinearConstraint(fun=Hc, lb=1, ub=1),
        NonlinearConstraint(fun=Lc, lb=1, ub=1),
    ),
)
assert res.success, res.message
print(res)

问题分析与修复建议 #

1. 无效值警告（RuntimeWarning） #

• 伽马函数简化：当z=1时，gamma(s[n1+i-1]+1)/gamma(s[n1+i-1]+1-z)可直接简化为s[n1+i-1]（因为gamma(a+1)/gamma(a)=a），替换后能避免伽马函数的数值不稳定问题。
• 幂运算异常处理：t**(s[n1+i]+i+1)中，若指数为负且t∈[0,1]，t=0时会出现0的负次幂，导致无穷大。可给t添加微小偏移（如t + 1e-8），或约束指数非负。

2. 迭代超限（AssertionError） #

• 放松约束范围：当前约束为固定值（如Sc必须等于43994），可行域极小甚至不存在。建议改为小范围约束，比如lb=43993, ub=43995，给优化器调整空间。
• 修正初始值维度：x0=np.ones(1 + 2*(n1+n2+n3+n4)-4)错误包含了t的维度，t是积分变量而非优化参数，初始值维度应与x一致（即2*(n1+n2+n3+n4)-4）。
• 更换优化器：默认的L-BFGS-B不适合复杂非线性约束问题，改用trust-constr优化器，同时调整迭代参数，比如增加maxiter、设置gtol。

3. 代码冗余清理 #

• 移除decimal精度设置：decimal.getcontext().prec=4会降低计算精度，导致数值不稳定，直接删除。
• 删除无用转换：S0=np.array(S)等未实际使用的转换代码，删除以减少不必要的计算。

修改后的关键代码片段 #

# 修正初始值维度，仅优化x变量
x_dim = 2*(n1+n2+n3+n4)-4
x0 = np.ones(x_dim)

# 放松约束为小范围
constraints = (
    NonlinearConstraint(fun=Sc, lb=43993, ub=43995),
    NonlinearConstraint(fun=Ic, lb=0.99, ub=1.01),
    NonlinearConstraint(fun=Hc, lb=0.99, ub=1.01),
    NonlinearConstraint(fun=Lc, lb=0.99, ub=1.01),
)

# 使用trust-constr优化器并调整参数
res = minimize(
    fun=objective,
    x0=x0,
    constraints=constraints,
    method='trust-constr',
    options={'maxiter': 1000, 'gtol': 1e-6}
)

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