如何才能最小化这个问题？

Question

我想最小化以下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


  

Answer 1

为了提高效率，将 I2 重写为纯 Numpy 表达式，以及 S、I、H、L。

t 和 x 有 13 个自由度。

import numpy as np
from scipy.optimize import minimize, NonlinearConstraint


def I2(x: np.ndarray) -> float:
    ab = np.array([21.1324865405187, 78.8675134594813])
    x0246 = x[0: 8: 2, np.newaxis]
    x1357 = x[1: 8: 2, np.newaxis]
    x8910 = x[8: 12, np.newaxis]
    x01, x23, x45, x67 = x0246 + x1357*(ab**(x8910 + 1) + 1)
    
    term1 = np.stack((
        x01*(+0.125*x23 + 0.0028) - 0.0043217,
        x23*(-0.0056*x45 - 0.125*x01 + 0.3778) - 0.029*x67,
        x23*(+0.0056*x45 - 0.35) + 0.0028*x45 - 0.535,
        -0.0250*x23 + 0.0318*x67,
    ))
    term2 = ab**x8910 * (x8910 + 1) * x1357
    return 50*((term1 + term2)**2).sum()



def solve() -> None:
    def objective(params: np.ndarray) -> float:
        x = params[1:]
        return I2(x)

    def S(params: np.ndarray) -> float:
        t, x = np.split(params, (1,))
        return t**(x[8] + 1)*x[1] + x[0] + x[1]
    def I(params: np.ndarray) -> float:
        t, x = np.split(params, (1,))
        return t**(x[9] + 1)*x[3] + x[2] + x[3]
    def H(params: np.ndarray) -> float:
        t, x = np.split(params, (1,))
        return t**(x[10] + 1)*x[5] + x[4] + x[5]
    def L(params: np.ndarray) -> float:
        t, x = np.split(params, (1,))
        return t**(x[11] + 1)*x[7] + x[6] + x[7]

    res = minimize(
        fun=objective,
        x0=np.ones(13),
        constraints=(
            NonlinearConstraint(fun=S, lb=43994, ub=43994),
            NonlinearConstraint(fun=I, lb=1, ub=1),
            NonlinearConstraint(fun=H, lb=0, ub=0),
            NonlinearConstraint(fun=L, lb=1, ub=1),
        ),
    )
    assert res.success, res.message
    print(res.x)


if __name__ == '__main__':
    solve()

[ 1.96620968e+05  2.12508470e+05 -1.68514470e+05  1.94073064e+05
 -1.94072064e+05  4.29304595e+04 -4.29304595e+04  3.93272854e+04
 -3.93262854e+04 -2.20148565e+06 -2.25335898e+06 -5.25227884e+04
 -4.09384848e+01]