My approach would be: first determine the time evolution of $\hat{x}(t)$ and $\hat{p}(t)$. For $\hat{x}$ you have
$$
\frac{d}{dt}\hat{x}_H(t) = i[H_H,\hat{x}_H(t)] = \frac{i}{2m} [\hat{p}_H(t)^2,\hat{x}_H(t)] = \frac{\hat{p_H(t)}}{m}
$$
and for $p$ you have (assuming $0\leq t \leq T$)
$$
\frac{d}{dt}\hat{p}_H(t) = i[H_H(t),\hat{p}_H(t)] = -m\omega_0^2 \hat{x}_H(t) + F_0\sin(\Omega t)
$$
These are coupled differential equations, which you can decouple by differentiating them once more with respect to time and performing a substitution. For instance,
$$
\frac{d^2}{dt^2} \hat{x}_H(t) = \frac{1}{m} \frac{d}{dt} \hat{p}_H = -\omega_0^2 \hat{x}_H(t)+\frac{F_0}{m}\sin(\Omega t)
$$
where I substituted $\frac{d}{dt} \hat{p}_H(t)$ by its equation of motion found earlier. You can also get an equation like this for $\hat{p}_H(t)(t)$, which I leave for you..
Now, these equations can be solved using your favorite method, provided you give them suitable boundary conditions. Note that you only need one boundary condition for $x$ and $p$ (which is $x_H(0)=\hat{x}_S$ and $p_H(0)=\hat{p}_S$ It will give you some expression for $\hat{x}_H(t)$ and $\hat{p}_H(t)$ in terms of $\hat{x}_S$ and $\hat{p}_S$. The Heisenberg Hamiltonian is then easily determined by substituting $\hat{x}_H(t)$ and $\hat{p}_H(t)$.
With that expression in hand you should be able to find $\langle H(t)\rangle$ (note that you should consider the cases where $t<0$ and $t>T$ separately).
EDIT:
The proof regarding my statement below: In the Schroedinger picture the Hamiltonian is
$$\hat H_S=\frac{\hat{p}_S^2}{2m}+\frac{1}{2}m\omega_0^2\hat{x}_S^2-\hat{x}_SF_0\sin(\Omega t)$$
and the Heisenberg picture is given by $H_H = U^\dagger(t) H_S U(t)$. So if you take for instance the first term you get:
$$U^\dagger(t) \frac{\hat{p}_S^2}{2m}U(t) =\frac{1}{2m} (U^\dagger(t) \hat{p}_S U^\dagger(t))(U(t)\hat{p}_SU(t)) =\frac{1}{2m} \hat{p}_H(t)^2 $$
You can do the same for the other terms. In the end you just effectively replace $p_S\rightarrow p_H(t)$ and the same for $x$.
First things first: the operator which corresponds to the energy is the Hamiltonian, typically written as $H$. So when you want to get the expectation value of the energy, you evaluate $\langle H\rangle$.
Now, there are multiple ways to do this. One way is to use the Schroedinger equation to get
$$\langle H\rangle = \left\langle i\hbar\frac{\partial}{\partial t}\right\rangle = \int\Psi^*(x,t) i\hbar\frac{\partial}{\partial t}\Psi(x,t)\,\mathrm{d}x\tag{1}$$
That calculation is completely general, i.e. it is valid in any situation for which the Schroedinger equation applies.
Another way to get $\langle H\rangle$ is to use the definition of the Hamiltonian operator, which in nonrelativistic QM is
$$H = \frac{p^2}{2m} + V(x,t) = -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x,t)$$
That gives you
$$\begin{align}\langle H\rangle &= \biggl\langle -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x,t)\biggr\rangle \\ &= \int\Psi^*(x,t)\biggl(-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x,t)\biggr)\Psi(x,t)\,\mathrm{d}x\tag{2}\end{align}$$
Either (1) or (2) works in general.
For a free particle only, the potential $V(x,t)$ is zero, and you get
$$\langle H\rangle = \int\Psi^*(x,t)\biggl(-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\biggr)\Psi(x,t)\,\mathrm{d}x$$
which is the expression you saw on that web page.
Best Answer
The key here is to notice that at any time $t$, one can expand the full time-dependent state $\Psi(x,t)$ at that time in terms of a basis $u_n$; $$ |\Psi(t)\rangle = \sum_n c_n(t)|u_n\rangle $$ The basic idea here is that the time-dependence of the state is being encoded in the expansion coefficients $c_n(t)$. Having done this, notice now that $$ \langle\Psi(t) |\hat A|\Psi(t)\rangle = \sum_{m,n}c_m^*(t)c_n(t)\langle u_m|\hat A|u_n\rangle = \sum_{m,n}c_m^*(t)A_{mn}c_n(t) $$
It is then possible to show (using the Schrodinger equation) that if the states $u_n$ are eigenstates of the Hamiltonian with corresponding eigenvalues $E_n$, then the coefficient $c_n(t)$ satisfy a simple differential equation in $t$ with a simple solution, and you should be able to put this all together to get what you're looking for.