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
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.