变量规模可变的Sympy方程组I2最小化方案求助

需要对以下Python代码中的I2进行最小化处理。尝试过scipy、sympy、gekko等库，但这些库要求显式定义变量，而变量数量会随代码中n的取值变化（如20、28、36个变量），显式定义大量变量并不合理。当前通过x = [sp.symbols('x%d' % i) for i in range(8*n-4)]生成变量，询问是否需要改用其他方法，或是对现有方法进行修改。

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

# create variables
n=2;z=1;m=2;a=0;b=100
x = [sp.symbols('x%d' % i) for i in range(8*n-4)]
#x=sp.symbols('x:20')
#t=sp.symbols('t ')
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]=factorial(i+j)/(2**j*factorial(j)*factorial(i-j))
PS=Matrix(1,n,x[0:n]);PI=Matrix(1,n,x[n:2*n]);PH=Matrix(1,n,x[2*n:3*n]);PL=Matrix(1,n,x[3*n:4*n])
TS=[[1]];TI=[[1]];TH=[[1]];TL=[[1]]
for i in range(n-1):
    TS.append([t**(x[4*n+i]+i+1)])
    TI.append([t**(x[5*n+i-1]+i+1)])
    TH.append([t**(x[6*n+i-2]+i+1)])
    TL.append([t**(x[7*n+i-3]+i+1)])
MTS=Matrix(TS);MTI=Matrix(TI);MTH=Matrix(TH);MTL=Matrix(TL)
S=PS*B*MTS;I=PI*B*MTI;H=PH*B*MTH;L=PL*B*MTL
#CONVERT SYMPY MATRICES TO NUMPY ONE
S0=np.array(S);I0=np.array(I);H0=np.array(H);L0=np.array(L)
GS=zeros(n,n);GI=zeros(n,n);GH=zeros(n,n);GL=zeros(n,n)
for i in range(n):
   if i==0:
      GS[i,i]=0
      GI[i,i]=0
      GH[i,i]=0
      GL[i,i]=0
   else:
      GS[i,i]=simplify(gamma(x[4*n+i-1]+i+1)/gamma(x[4*n+i-1]+i+1-z))
      GI[i,i]=simplify(gamma(x[5*n+i-2]+i+1)/gamma(x[5*n+i-2]+i+1-z))
      GH[i,i]=simplify(gamma(x[6*n+i-3]+i+1)/gamma(x[6*n+i-3]+i+1-z))
      GL[i,i]=simplify(gamma(x[7*n+i-4]+i+1)/gamma(x[7*n+i-4]+i+1-z))
DS=simplify(t**(-z)*PS*B*GS*MTS)
DI=simplify(t**(-z)*PI*B*GI*MTI)
DH=simplify(t**(-z)*PH*B*GH*MTH)
DL=simplify(t**(-z)*PL*B*GL*MTL)
RS=DS[0,0]-(.0043217-.125*S[0,0]*I[0,0]-(.002+.0008)*S[0,0])
RI=DI[0,0]-(.125*S[0,0]*I[0,0]+.0056*H[0,0]*I[0,0]+.029*L[0,0]-(.002+.0008+.025+.35)*I[0,0])
RH=DH[0,0]-(.535-.0056*H[0,0]*I[0,0]+.35*I[0,0]-(.002+.0008)*H[0,0])
RL=DL[0,0]-(.025*I[0,0]-(.002+.0008+.029)*L[0,0])
#f = lambdify(t, R)
def R(s):
   return (RS**2+RI**2+RH**2+RL**2).subs(t,s)
r, w = roots_legendre(m)
w=zeros(1,m)
I1=0
for i in range (m):
   w[i]=R((b-a)/2*r[i]+(b+a)/2)
   I1=I1+w[i]
I2=((b-a)/2)*I1

不需要放弃当前的变量生成方式，可通过以下两种适配方案实现最小化：

1. 适配Scipy优化器：符号转数值函数 #

利用Sympy的lambdify工具，将符号表达式I2转换为接受numpy数组输入的数值函数，直接适配Scipy的优化接口，无需显式定义单个变量：

# 先确保I2是仅依赖x的符号表达式（t已被积分消除）
# 将I2转为支持数组输入的数值函数
I2_func = sp.lambdify(x, I2, 'numpy')

# 定义优化目标函数，接受一维numpy数组
def objective(x_arr):
    return I2_func(*x_arr)

# 生成初始值数组，长度与变量列表一致
x0 = np.ones(len(x))

# 调用Scipy优化器求解
from scipy.optimize import minimize
result = minimize(objective, x0, method='L-BFGS-B')

print("最优变量值：", result.x)
print("最小I2值：", result.fun)

2. 适配Gekko：动态生成变量 #

Gekko支持动态批量创建变量，无需逐个显式定义，只需将原Sympy逻辑转换为Gekko支持的运算语法即可：

from gekko import GEKKO

# 初始化本地Gekko模型
m = GEKKO(remote=False)
# 动态生成变量列表，数量随n自动变化
x = [m.Var(value=1) for _ in range(8*n-4)]

# 以下为原Sympy逻辑的Gekko转换示例（需替换对应运算）
# 1. 构建B矩阵
B = m.Array(m.Const, (n,n))
for i in range(n):
    for j in range(n):
        if i < j:
            B[i,j] = 0
        else:
            # Gekko中阶乘可用m.factorial
            B[i,j] = m.factorial(i+j)/(2**j * m.factorial(j) * m.factorial(i-j))

# 2. 构建PS、PI等矩阵（Gekko矩阵操作）
PS = m.Array(m.Var, (1,n))
for i in range(n):
    PS[0,i] = x[i]
# 其余PI、PH、PL同理...

# 3. 处理gamma函数（Gekko直接支持m.gamma）
# 其余逻辑按Gekko语法转换后，定义目标函数
m.Obj(I2)

# 求解
m.solve(disp=True)

3. 可选优化：减少变量规模 #

观察代码中变量的使用模式，例如TS、TI中的指数变量是否存在可合并的约束关系，可尝试归类变量，减少总变量数，降低优化问题的复杂度。

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