How to minimize a numerical integration in python?

Question

I need some help to minimize a numerical integration. It's about a classical problem in physics (hydrogen atom). It can be solved analytically but I need to solve it numerically in Python.

We have an electron's wave function $\psi(r,a)=e^{-r/a}$ and a Hamiltonian (in Ry units) $H=-\nabla^2-\frac{2}{r}$. The expected value of Energy is given by

$$\langle E(a)\rangle =\frac{\iiint \psi^{*}(r,a) H \psi(r,a) r^2 \sin(\theta) \mathrm dr \mathrm d\theta \mathrm d\phi}{\iiint |\psi(r,a)|^2 r^2 \sin(\theta) \mathrm dr \mathrm d\theta \mathrm d\phi}$$.

The integral it's over all space ($0\leq r \leq \infty$ ; $0\leq \theta \leq \pi$; $0\leq \phi \leq 2\pi$ )

Doing $\langle E \rangle$ numerically, how to find (also numerically) the value of $a$ which minimizes the value of the integral?

Edit 1- I know the value of '$a$' that minimizes the value of $<E(a)>$.

$a=1$ and $<E(1)>=-1$. But somehow I need to 'discovery' the value of $'a'$ that minimizes $<E(a)>$ numerically.

The code:

import numpy as np
from scipy.integrate import tplquad
import sympy as sp

# deffining the symbols
r, theta, phi, a = sp.symbols('r theta phi a')

# Deffining the function
f = sp.exp(-r/a)

# laplacian
laplacian = sp.diff(r**2 * sp.diff(f, r), r) / r**2 + \
            sp.diff(sp.sin(theta) * sp.diff(f, theta), theta) / 
(r**2 * sp.sin(theta)) + \
            sp.diff(f, phi, phi) / (r**2 * sp.sin(theta)**2)

# changing the functions into a lambda functions
laplacian_numeric = sp.lambdify((r, theta, phi, a), laplacian, 
'numpy')

f_numeric = sp.lambdify((r, theta, phi, a), f, 'numpy')

# deffining the kernel of integrations
def integrand(r, theta, phi, a):
    jac = r**2 * sp.sin(theta)
    return f_numeric(r, theta, phi, a) * laplacian_numeric(r, 
theta, phi, a)*jac

def integrand2(r,theta,phi,a):
    jac=r**2 * sp.sin(theta)
    return 2*jac*f_numeric(r, theta, phi, a)**2 /r

def integrand3(r,theta,phi,a):
    jac=r**2 * sp.sin(theta)
    return jac*f_numeric(r, theta, phi, a)**2

# integrations limits
r_limits = (0, np.inf)
theta_limits = (0, np.pi)
phi_limits = (0, 2*np.pi)

# value of the parameter a
a_value = 1

# numerical integration
h1, _ = tplquad(integrand, *phi_limits, *theta_limits,*r_limits, args=(a_value,))

h2,_ = tplquad(integrand2,*phi_limits, *theta_limits,*r_limits, args=(a_value,))

norm,_ = tplquad(integrand3,*phi_limits, *theta_limits,*r_limits, args=(a_value,))

print('result of numerical integration', (-h1-h2)/norm)