pythonoptimizationminimize

Minimize sympy system of equations of variable size


I want to minimize I2 in the following Python code. I tried to use some libraries like scipy, sympy, gekko ,... but all of them need to introduce variables explicitly while the number of variables changes with changing n in the code. On the other writing too many variables (for example 20, 28 ,36, ... variables) explicitly is illogical. In this code x = [sp.symbols('x%d' % i) for i in range(8*n-4)] is my variables. Should another method be used? Should this method be modified? I need your help.

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


  

Solution

  • For efficiency at n=2, rewrite I2 as a pure-Numpy expression, as well as S, I, H, L.

    With t and x there are 13 degrees of freedom.

    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]
    

    Fast and effective. Compare what you'd have to do to support a Sympy expression of arbitrary size; there are several different methods, one being

    t = next(s for s in S.free_symbols if s.name == 't')
    x = [
        next(s for s in I2.free_symbols if s.name == f'x{i}')
        for i in range(8*n - 4)
    ]
    args = [t, *x]
    I2 = sympy.lambdify(args, I2)
    S = sympy.lambdify(args, S.flat()[0])
    I = sympy.lambdify(args, I.flat()[0])
    H = sympy.lambdify(args, H.flat()[0])
    L = sympy.lambdify(args, L.flat()[0])
    
    def objective(params: np.ndarray) -> float:
        return I2(*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 + 8*n-4),
        constraints=(
            NonlinearConstraint(fun=Sc, lb=43994, ub=43994),
            NonlinearConstraint(fun=Ic, lb=1, ub=1),
            NonlinearConstraint(fun=Hc, lb=0, ub=0),
            NonlinearConstraint(fun=Lc, lb=1, ub=1),
        ),
    )
    assert res.success, res.message
    

    Much slower.