How to Put an indefinite Integral in an objective function Python

Question

I have this objective function that I want to minimize in regards to 't', and in this function, there is an integral from 0 to t that has to be considered for the optimization. However, I have no Idea how to include this integral to my objective function. I've used simpy and scipy library for that but none of them seems to work. here is my code:

Cf= 100;
Cp=50;
  Eta=72.66585511711865;
  Beta=1.18609324
def R(t):
  return math.exp(-t/Eta)**Beta
def f(t):
  return (Beta/(Eta**Beta))*t**(Beta-1)*(math.exp(-t/Eta)**Beta)  
def M(t):
  return t*f(t)
import sympy as sp 
from scipy import *
import scipy.integrate
t = sp.Symbol('t') 
Mt=t*(Beta/(Eta**Beta))*t**(Beta-1)*(sp.exp(-t/Eta)**Beta);
intM=sp.integrate(Mt,t)
print(intM)
def C(t):
  return (Cp*R(t)+Cf*(1-R(t)))/(t*R(t)+ intM)
from scipy.optimize import minimize_scalar
res = minimize_scalar(C)
res

this gives me error " TypeError: cannot determine truth value of Relational" . in scipy.integrate.quad, I can't define my upper limit as 't'. And on the other hand, it won't work when I try to put the sp.integrate output to the objective function. I would be happy if u could help. thanks.

Answer 1

First up, you can't have an indefinite integral within an optimization routine, as that's defined up to a constant. Not a function, but a class of functions.

But it looks like yours is a definite integral - you integrate it over some dummy variable, say, t', from t'=0 to t'=t. So that's perfectly fine, and you just have to do it.

First, let's do your definite integral:

intM = sp.integrate(Mt, (t, 0, t))
# second arg is variable of integration, lower bound, upper bound
# we can skip defining tprime, a mathematician would hate you for this but meh

This still depends on t, so now let's turn it into a numerical function:

intM_func = sp.lambdify(t, intM)
# note how you can do intM_func(3) and get a number

Then stick intM_func(t) rather than the symbolic intM into your C minimizable function and you should be good to go :)

---

The full version that works for me:

Cf = 100
Cp = 50
Eta = 72.66585511711865
Beta = 1.18609324


def R(t):
    return math.exp(-t / Eta) ** Beta


def f(t):
    return (Beta / (Eta**Beta)) * t ** (Beta - 1) * (math.exp(-t / Eta) ** Beta)


def M(t):
    return t * f(t)


import sympy as sp
from scipy import *
import scipy.integrate

t = sp.Symbol("t")
Mt = t * (Beta / (Eta**Beta)) * t ** (Beta - 1) * (sp.exp(-t / Eta) ** Beta)
intM = sp.integrate(Mt, (t, 0, t))
intMfunc = sp.lambdify(t, intM)
print(intM)


import numpy as np
def C(t):
    if t == 0:   # dodge a division-by-zero bullet
        return np.inf
    return (Cp * R(t) + Cf * (1 - R(t))) / (t * R(t) + intMfunc(t))


from scipy.optimize import minimize_scalar

res = minimize_scalar(C)
print(res)

And I get:

     fun: 1.5407809938973502
 message: '\nOptimization terminated successfully;\nThe returned value satisfies the termination criteria\n(using xtol = 1.48e-08 )'
    nfev: 61
     nit: 28
 success: True
       x: 2964.8614509894874

as an output.

Answer 2

I tempted to 'close' this because you did not give the full error - with traceback. But it appears you got an answer, so I'll just try to focus the error as you gave it.

error " TypeError: cannot determine truth value of Relational" . in scipy.integrate.quad

This is similar to a common "ambiguity error" with numpy arrays. Something is trying to do a simple True/False action, such as a if. The scipy function might be testing whether the range values are in correct order. Or maybe the error occurs in your func. Without traceback we are left to guessing, or deducing from previous knowledge. So you are lucky that someone was already familiar with these issues. Just don't count on that!

To illustrate:

In an isympy session, consider a symbol:

In [35]: x

 

And a relational:

In [36]: x>0

In [37]: type(_)
Out[37]: sympy.core.relational.StrictGreaterThan

Putting that relational in an `if:

In [38]: if x>0: print('yes')
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [38], in <cell line: 1>()
----> 1 if x>0: print('yes')

File ~\anaconda3\lib\site-packages\sympy\core\relational.py:511, in Relational.__bool__(self)
    510 def __bool__(self):
--> 511     raise TypeError("cannot determine truth value of Relational")

TypeError: cannot determine truth value of Relational

sympy can't return a simple True/False from a Relational expression.

Using sympy and scipy together requires considerable knowledge of how both work. It's all too easy to using sympy expressions in a numeric function in ways that don't work.

Here's the missing traceback:

In [59]: res = minimize_scalar(C)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [59], in <cell line: 1>()
----> 1 res = minimize_scalar(C)

File ~\anaconda3\lib\site-packages\scipy\optimize\_minimize.py:794, in minimize_scalar(fun, bracket, bounds, args, method, tol, options)
    792     return method(fun, args=args, bracket=bracket, bounds=bounds, **options)
    793 elif meth == 'brent':
--> 794     return _minimize_scalar_brent(fun, bracket, args, **options)
    795 elif meth == 'bounded':
    796     if bounds is None:

File ~\anaconda3\lib\site-packages\scipy\optimize\optimize.py:2396, in _minimize_scalar_brent(func, brack, args, xtol, maxiter, **unknown_options)
   2393 brent = Brent(func=func, args=args, tol=tol,
   2394               full_output=True, maxiter=maxiter)
   2395 brent.set_bracket(brack)
-> 2396 brent.optimize()
   2397 x, fval, nit, nfev = brent.get_result(full_output=True)
   2399 success = nit < maxiter and not (np.isnan(x) or np.isnan(fval))

File ~\anaconda3\lib\site-packages\scipy\optimize\optimize.py:2180, in Brent.optimize(self)
   2177 def optimize(self):
   2178     # set up for optimization
   2179     func = self.func
-> 2180     xa, xb, xc, fa, fb, fc, funcalls = self.get_bracket_info()
   2181     _mintol = self._mintol
   2182     _cg = self._cg

File ~\anaconda3\lib\site-packages\scipy\optimize\optimize.py:2154, in Brent.get_bracket_info(self)
   2151 ### BEGIN core bracket_info code ###
   2152 ### carefully DOCUMENT any CHANGES in core ##
   2153 if brack is None:
-> 2154     xa, xb, xc, fa, fb, fc, funcalls = bracket(func, args=args)
   2155 elif len(brack) == 2:
   2156     xa, xb, xc, fa, fb, fc, funcalls = bracket(func, xa=brack[0],
   2157                                                xb=brack[1], args=args)

File ~\anaconda3\lib\site-packages\scipy\optimize\optimize.py:2606, in bracket(func, xa, xb, args, grow_limit, maxiter)
   2604 fa = func(*(xa,) + args)
   2605 fb = func(*(xb,) + args)
-> 2606 if (fa < fb):                      # Switch so fa > fb
   2607     xa, xb = xb, xa
   2608     fa, fb = fb, fa

File ~\anaconda3\lib\site-packages\sympy\core\relational.py:511, in Relational.__bool__(self)
    510 def __bool__(self):
--> 511     raise TypeError("cannot determine truth value of Relational")

TypeError: cannot determine truth value of Relational

So as I suspected, scipy is trying to check ranges;

if (fa < fb):                      # Switch so fa > fb

C is given a simple scalar argument produces a sympy expression - with an unevaluated integral.

In [57]: C(1)
Out[57]: 

That intM that you forgot to show us:

In [60]: intM
Out[60]: