Minimize sympy system of equations of variable size

Question

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

Please clarify the actual problem statement. What is the definition of the function that you want to minnimize. Some mathematical formula might be beneficial — FlyingTeller, Aug 15 '23 at 12:05
I want to minimize "I2" which was shown at the end of this code. "I2" is the answer of an integral — hojat saeidi, Aug 15 '23 at 12:11
`I2` is a function of `x0` through `x11`. When you minimize it with respect to `x`, are there bounds on `x`? Are there any other constraints? — Reinderien, Aug 15 '23 at 13:09
there is no bound for `x` and there are 4 constraints for `I2` : `S(0)-43994=0` , `I(0)-1=0` , `H(0)=0` and `L(0)-1=0` . — hojat saeidi, Aug 15 '23 at 13:39
I2 is not defined in terms of S, H and L; it's only defined in terms of x. If there's a secondary relationship between those variables you should explain that. — Reinderien, Aug 16 '23 at 01:56
In fact, all of the variables that you referenced in your constraints are functions of `t`. I have no option but to assume that that's a free, continuous variable. — Reinderien, Aug 16 '23 at 02:01

Reinderien · Accepted Answer · 2023-08-16T06:18:53.223

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.

Thanks for your answer @Reinderien but I don't understand why do you define `ab = np.array([21.1324865405187, 78.8675134594813])`. from where do you get the numbers `21.1324865405187, 78.8675134594813` ?In `term1` the varable`x01` exactly refers to `S` in my code but in your code there is no relationship between `x01` and `S`; you defined them separately. I can't understand what do you define ` x01, x23, x45, x67 = x0246 + x1357*(ab**(x8910 + 1) + 1)` because I can't find any relation between this and my code. I would be grateful if you could explain — hojat saeidi, Aug 16 '23 at 08:58
`ab` is a consequence of your expression for `I2` and is converted manually for demonstration (refer to the output of `print(I2)` to see where those terms appear in the expanded form). — Reinderien, Aug 16 '23 at 15:17
_in your code there is no relationship between x01 and S_ - that's not strictly correct. There is a direct mathematical relationship, but the expression is a manual conversion for the case of n=2. For the general case refer to the second block of code. — Reinderien, Aug 16 '23 at 15:21
_I can't find any relation between [the expression for x01] and my code_ - it's mathematically equivalent. Do the integral analytically by hand or examine the structure of `I2` and you will find as much — Reinderien, Aug 16 '23 at 15:26
For `n=2` we have 12 unknown variables(8n-4) but in your output code there are 14. why? — hojat saeidi, Aug 16 '23 at 19:34
There aren't 14, there are 13. `t` is unknown, like I said in the comments on the question. — Reinderien, Aug 16 '23 at 21:22
Yes, I made a mistake in counting the variables but `t` is a variable that we integrate from `I2` with respect to it (in interval [a,b] ) and must be eliminated, so we must have 12 variables. With this assumption, how should your code be modified? — hojat saeidi, Aug 17 '23 at 12:28
If you want to constrain S, I, H, L in the same parameter space they also need to be integrated. Have you shown integral forms of these variables in your question? — Reinderien, Aug 17 '23 at 13:05
The integral of S, I, H, L exist implicitly in the formula `RS**2+RI**2+RH**2+RL**2` . in fact i don't need the integral of these function, i need the integral of remainder that shown in form `I2` — hojat saeidi, Aug 17 '23 at 17:28

Minimize sympy system of equations of variable size

1 Answers1