Evaluate $\int_0^{\pi/4}{(4\cot x\ln\sec x-x)\ln^2\tan xdx}$
Here is my attempt at a solution. I did not get your very nice answer of $\frac{5\pi^4}{2304}$, but if yours and my answers are equivalent, it will mean we have found a value for $\operatorname{Re} \operatorname{Li}_4 (1 + i)$ (is a closed-form value known for this quantity?).
Set $$I = \int_0^{\frac{\pi}{4}} \ln^2 \tan x (4 \cot x \ln \sec x - x) \, dx.$$ Enforcing a substitution of $x \mapsto \arctan x$ leads to $$I = 2 \int_0^1 \frac{\ln^2 x}{x} \frac{\ln (1 + x^2)}{1 + x^2} \, dx - \int_0^1 \frac{\ln^2 x \arctan x}{1 + x^2} \, dx = 2I_1 - I_2.$$
First integral $I_1$
Making use of the following generating function for the harmonic numbers, namely $$\frac{\ln (1 + x^2)}{1 + x^2} = -\sum_{n = 1}^\infty (-1)^n H_n x^{2x}.$$ we have \begin{align} I_1 &= -\sum_{n = 1}^\infty (-1)^n H_n \int_0^1 x^{2n - 1} \ln^2 x \, dx\\ &= -\sum_{n = 1}^\infty (-1)^n H_n \frac{d^2}{ds^2} \left [\int_0^1 x^{2n + s -1} \, dx \right ]_{s = 0}\\ &= -\sum_{n = 1}^\infty (-1)^n H_n \frac{d^2}{ds^2} \left [\frac{1}{2n + s} \right ]_{s = 0}\\ &= -\frac{1}{4} \sum_{n = 1}^\infty \frac{(-1)^n H_n}{n^3}. \end{align}
Second integral $I_2$
Taking the Cauchy product between the Maclaurin series expansions for $\arctan x$ and $\frac{1}{1 + x^2}$ one finds $$\frac{\arctan x}{1 + x^2} = \sum_{n = 0}^\infty (-1)^n \left (H_{2n + 1} - \frac{1}{2} H_n \right ) x^{2n + 1}.$$ Thus \begin{align} I_2 &= \sum_{n = 0}^\infty (-1)^n \left (H_{2n + 1} - \frac{1}{2} H_n \right ) \int_0^1 x^{2n + 1} \ln^2 x \, dx\\ &= \sum_{n = 0}^\infty (-1)^n \left (H_{2n + 1} - \frac{1}{2} H_n \right ) \frac{d^2}{ds^2} \left [\int_0^1 x^{2n + s + 1} \, dx \right ]_{s = 0}\\ &= \sum_{n = 0}^\infty (-1)^n \left (H_{2n + 1} - \frac{1}{2} H_n \right ) \frac{d^2}{ds^2} \left [\frac{1}{2n + s + 2} \right ]_{s = 0}\\ &= \frac{1}{4} \underbrace{\sum_{n = 0}^\infty \frac{(-1)^n}{(n + 1)^3} \left (H_{2n + 1} - \frac{1}{2} H_n \right )}_{n \, \mapsto \, n - 1}\\ &= \frac{1}{4} \sum_{n = 1}^\infty \frac{(-1)^{n - 1}}{n^3} \left (H_{2n - 1} - \frac{1}{2} H_{n - 1} \right )\\ &= -\frac{1}{4} \sum_{n = 1}^\infty \frac{(-1)^n}{n^3} \left (H_{2n} - \frac{1}{2n} \right ) + \frac{1}{8} \sum_{n = 1}^\infty \frac{(-1)^n}{n^3} \left (H_n - \frac{1}{n} \right )\\ &= -\frac{1}{4} \sum_{n = 1}^\infty \frac{(-1)^n H_{2n}}{n^3} + \frac{1}{8} \sum_{n = 1}^\infty \frac{(-1)^n H_n}{n^3} \end{align}
Main integral $I$
So for the main integral $I$ we have $$I = -\frac{5}{8} \sum_{n = 1}^\infty \frac{(-1)^n H_n}{n^3} + \frac{1}{4} \sum_{n = 1}^\infty \frac{(-1)^n H_{2n}}{n^3}.$$
Dealing with these two Euler sums, their values can be found from the following generating function \begin{align} \sum^\infty_{n=1}\frac{H_n}{n^3}x^n &=2{\rm Li}_4(x)+{\rm Li}_4\left(\tfrac{x}{x-1}\right)-{\rm Li}_4(1-x)-{\rm Li}_3(x)\ln(1-z)-\frac{1}{2}{\rm Li}_2^2\left(\tfrac{x}{x-1}\right)\\ &+\frac{1}{2}{\rm Li}_2(x)\ln^2(1-x)+\frac{1}{2}{\rm Li}_2^2(x)+\frac{1}{6}\ln^4(1-x)-\frac{1}{6}\ln{x}\ln^3(1-x)\\ &+\frac{\pi^2}{12}\ln^2(1-x)+\zeta(3)\ln(1-x)+\frac{\pi^4}{90},\tag1 \end{align} which is proved in this answer here.
Setting $x = -1$ in (1) gives \begin{align} \sum^\infty_{n=1}\frac{(-1)^nH_n}{n^3}=2{\rm Li}_4\left(\tfrac{1}{2}\right)-\frac{11\pi^4}{360}+\frac{7}{4}\zeta(3)\ln{2}-\frac{\pi^2}{12}\ln^2{2}+\frac{1}{12}\ln^4{2}, \end{align} while setting $x = i$ in (1) gives \begin{align} \frac{1}{4} \sum_{n = 1}^\infty \frac{(-1)^n H_{2n}}{n^3} &= 2 \sum_{n = 1}^\infty \frac{(-1)^n H_{2n}}{(2n)^3}\\ &= 2 \operatorname{Re} \sum_{n = 1}^\infty \frac{H_n}{n^3} i^n\\ &= -4 \operatorname{Re} \operatorname{Li}_4(1 + i) + \frac{29 \pi^4}{1152} + \frac{35}{32} \zeta (3) \ln 2 + \frac{\pi^2}{32} \ln^2 2. \end{align}
Substituting these two values for the Euler sums back into the expression for the integral $I$ gives a final answer of
$$I = -\frac{5}{4} \operatorname{Li}_4 \left (\frac{1}{2} \right ) - 4 \operatorname{Re} \operatorname{Li}_4 (1 + i) + \frac{17}{384} \pi^4 + \frac{\pi^2}{12} \ln^2 2 - \frac{5}{96} \ln^4 2.$$
So, this leads one to the following conjecture. Does?
$$\operatorname{Re} \operatorname{Li}_4 (1 + i) = -\frac{5}{16} \operatorname{Li}_4 \left (\frac{1}{2} \right ) + \frac{97}{9216} \pi^4 + \frac{\pi^2}{48} \ln^2 2 - \frac{5}{384} \ln^4 2$$
Update
The conjecture is true! A proof of this can be found here. So one does indeed have $$\int_0^{\frac{\pi}{4}} \ln^2 \tan x (4 \cot x \ln \sec x - x) \, dx = \frac{5 \pi^4}{2304}.$$ It would of course be nice to find a simple (simpler?) solution to this integral that, unlike my solution, does not rely on knowing the value for $\operatorname{Re} \operatorname{Li}_4 (1 + i)$.
$\DeclareMathOperator{\arctanh}{arctanh}$ We may write \begin{align*} I =&\int_0^{\frac \pi 4} \Big(4\cot x\log (\sec x) -x \Big)\log^2 (\tan x) \ dx\\ =& 4\underbrace{\int_0^{\frac \pi 4} \cot x\log (\sec x) \log^2 (\tan x) \ dx}_{\tan x \mapsto x} -\underbrace{\int_0^{\frac \pi 4} x \log^2 (\tan x) \ dx}_{=:I_2}\\ =& 2\underbrace{\int_0^1 \frac{\log^2 x\log(1+x^2)}{x(1+x^2)} dx}_{x^2 \to x} - I_2\\ =&\underbrace{\frac 1 4 \int_0^1 \frac{\log^2 x \log(1+x)}{x(1+x)}dx}_{=:I_1} - I_2\\ =&I_1 - I_2. \end{align*} Using the Maclaurin series of $\frac{\log(1+x)}{1+x} = -\sum_{k=1}^\infty (-1)^{k}H_k x^k,$ we get \begin{align*}I_1 =&- \frac 1 4 \sum_{k=1}^\infty (-1)^{k}H_k \int_0^1 x^{k-1} \log^2 x\ dx\\ =&\boxed{-\frac 1 2 \sum_{k=1}^\infty \frac{(-1)^{k}H_k}{k^3}} \end{align*} as can be found in @omegadot's answer. (In fact, it holds that $$ \sum_{k=1}^\infty \frac{(-1)^{k}H_k}{k^3}=2\text{Li}_4(1/2) + \frac {7\ln 2 \zeta(3)} 4 - \frac {\pi^2 \ln^2 2 }{12} -\frac {11\pi^4}{360} + \frac {\ln^4 2}{12}, $$ though we don't really need this fact for our purpose now.)
Next, we derive that: \begin{align*} \boxed{I_2 = \frac 1 8 \int_0^1 \frac{\log x \log^2\left(\frac {1+x}{1-x}\right)}{x}dx +\frac{\pi^4}{256}.} \end{align*} Then, using the Maclaurin series of $$\displaystyle\frac 1 4 \log^2\left(\frac {1+x}{1-x}\right)=\arctanh^2 x =\sum_{k= 1}^\infty \frac {H_{2k}-\frac 1 2 H_k}{k}x^{2k} $$ we get: \begin{align*} \frac 1 8 \int_0^1 \frac{\log x \log^2\left(\frac {1+x}{1-x}\right)}{x}dx =&\sum_{k=1}^\infty \frac{H_{2k}-\frac 1 2 H_k}{2k}\int_0^1 x^{2k-1} \log x\ dx\\ =- &\sum_{k=1}^\infty \frac{H_{2k}-\frac 12 H_k}{8k^3} \\ =&-\frac 1 2\sum_{k=1}^\infty\frac{(1+(-1)^k)H_k}{k^3} +\frac 1 {16} \sum_{k=1}^\infty\frac{H_k}{k^3}\\ =&-\frac 1 2\sum_{k=1}^\infty\frac{(-1)^kH_k}{k^3} - \frac 7 {16} \sum_{k=1}^\infty\frac{H_k}{k^3}\\ =& I_1 - \frac 7 {16} \sum_{k=1}^\infty\frac{H_k}{k^3}. \end{align*} This implies that \begin{align*} I_2 =& I_1 - \frac 7 {16} \sum_{k=1}^\infty\frac{H_k}{k^3} +\frac{\pi^4}{256}\\ \end{align*} hence $$ \boxed{I = I_1 - I_2 = \frac 7{16}\sum_{k=1}^\infty\frac{H_k}{k^3}-\frac {\pi^4}{256} = \frac 7 {16}\frac {\pi^4} {72} - \frac {\pi^4}{256}= \frac{5\pi^4}{2304}.} $$ Here, Euler's formula $$ 2\sum_{k=1}^\infty\frac{H_k}{k^n} = (n+2)\zeta(n+1) - \sum_{k=1}^{n-2}\zeta(k+1)\zeta(n-k),\qquad n\ge 2$$ can be used to evaluate $\sum_{k=1}^\infty\frac{H_k}{k^3} = \frac{\pi^4}{72}$.
Evaluation of $I_2$: We note by the symmetry that \begin{align*} \int_0^{\frac \pi 4} \log^2(\tan(x)) dx= \frac 1 2\int_0^{\frac \pi 2} \log^2\left(\tan\left(\frac x 2\right)\right) dx= \frac 1 4 \int_{0}^{\pi} \log^2\left(\tan\left(\frac x 2\right)\right) dx = \frac{\pi^3}{16}. \end{align*} This follows from $\displaystyle \log\left(\tan\left(\frac x 2\right)\right) = \sum_{k=1}^\infty \frac{(-1)^k -1}{k} \cos(kx)$ and Parseval's identity. Thus it follows \begin{align*} I_2 =& \underbrace{\int_0^{\frac\pi 4} \left(x-\frac \pi 4\right) \log^2(\tan(x)) dx }_{\frac \pi 4 -x\mapsto x}+ \frac{\pi^4}{64}\\ =& -\underbrace{\int_0^{\frac \pi 4} x \log^2\left(\frac {1+\tan x}{1-\tan x}\right) dx}_{\tan x = y} + \frac{\pi^4}{64} \\ =& -\underbrace{4\int_0^1 \frac{\arctan y}{1+y^2}\arctanh^2 y\ dy}_{=:J}+ \frac{\pi^4}{64}\\=&-J +\frac{\pi^4}{64}. \end{align*}
Now let us define $$ f(z) = \frac{\arctan z}{1+z^2}\left( \arctanh^2 z +\frac {\pi^2} {16}\right) $$ so that $f$ is analytic on $|z|<1$ and continuous on $|z|\le 1$. We will use the symmetric nature of $f$ to evaluate the $J$.
We first notice that $$ \int_{0}^1 f(z)dz =\frac 1 4 J +\frac{\pi^2}{16}\int_0^1\frac{\arctan z}{1+z^2}dz =\frac 1 4 J + \frac{\pi^4}{512} . $$ On the other hand, we have that \begin{align*} \int_{0}^i f(z)dz=&i\int_0^1 f(ix)dx \\ =&i\int_0^1 \frac{\arctan(ix)}{1-x^2}\left(\arctanh^2(ix) + \frac{\pi^2}{16}\right) dx\\ =& \int_0^1 \frac{\arctanh x}{1-x^2} \left(\arctan^2 x - \frac{\pi^2}{16}\right) dx\\ \underset{\text{IBP}}{=}&-\int_0^1 \arctanh^2 x\ \frac{\arctan x}{1+x^2} dx\\ =&-\frac 1 4 J \end{align*} where we have used the relations $\arctan(ix) = i\arctanh x$ and $\arctanh(ix) = i \arctan x$. Combining these, we know that \begin{align*} J =& \frac 12 J +\frac 1 2 J\\ =&2\left(\int_{0}^1f(z)dz -\frac{\pi^4}{512}\right) -2 \int_{0}^i f(z)dz\\=&-\underbrace{2 \int_1^i f(z)dz}_{=:K} -\frac{\pi^4}{256}\\ =&-K-\frac{\pi^4}{256}. \end{align*}
For $K$, we parametrize $[1,i]$ by $z=e^{i\theta}, 0\le \theta \le \frac \pi 2$. Then exploiting the fact that $$\arctan(e^{i\theta}) = -\frac i 2 \Big[\log\big(\cot\big(\scriptsize{\frac \theta 2+\frac\pi 4} \normalsize\big)\big)+ \frac {\pi i} 2\Big]$$ $$ \arctanh(e^{i\theta}) = \frac 1 2 \Big[\log\big(\cot\big(\scriptsize \frac \theta 2 \normalsize\big)\big) +\frac{\pi i}2\Big] $$ we can see that (by Cauchy's integral theorem) \begin{align*} K=&2 \int_0^{\frac \pi 2} f(e^{i\theta}) ie^{i\theta} d\theta \\ =&2i\int_0^{\frac \pi 2} \frac{\arctan(e^{i\theta})}{1+e^{i2\theta}} \left(\arctanh^2(e^{i\theta}) + \frac{\pi^2}{16}\right)e^{i\theta}d\theta \\ =&\frac 1 {8} \int_0^{\frac\pi 2} \frac{\log\big(\cot\big(\frac \theta 2 +\frac \pi 4\big)\big)+\frac{\pi i}2}{\cos \theta}\Big[\log^2\big(\cot\big(\scriptsize \frac \theta 2 \normalsize\big)\big)+\pi i \log\big(\cot\big(\scriptsize \frac \theta 2 \normalsize\big)\big)\Big]d\theta. \end{align*} Since $K$ is real, we can take the real part of the integral and perform the change of variable $\tan(\frac \theta 2) = y$ to get \begin{align*} K =&\frac 1 {8}\int_0^{\frac\pi 2} \frac{\log^2\big(\cot\big(\frac \theta 2\big)\big)\log\big(\cot\big(\frac \theta 2 +\frac \pi 4\big)\big) -\frac{\pi^2}2 \log\big(\cot\big(\frac\theta 2\big)\big)}{\cos \theta}d\theta\\ =&\frac 1 {4} \int_0^1 \frac{1}{1-y^2}\left(\log^2 y\log\left(\frac {1-y}{1+y}\right) +\frac{\pi^2}2 \log y\right)dy \\ =&\frac 1 {4} \underbrace{\int_0^1 \frac{\log^2 y\log\left(\frac {1-y}{1+y}\right)}{1-y^2}dy}_{\frac{1-y}{1+y} = x} +\frac{\pi^2}{8} \int_0^1 \frac{\log y}{1-y^2}dy\\ =&\frac 1 {8}\int_0^1 \frac{\log x\log^2\left(\frac {1+x}{1-x}\right)}{x}dx +\frac{\pi^2}{8} \sum_{n=0}^\infty \int_0^1 y^{2n}\log y\ dy \\ =&\frac 1 {8}\int_0^1 \frac{\log x\log^2\left(\frac {1+x}{1-x}\right)}{x}dx -\frac{\pi^2}{8} \underbrace{\sum_{n=0}^\infty\frac 1{(2n+1)^2}}_{=\frac 3 4 \zeta(2) = \frac{\pi^2}{8}} \\ =&\frac 1 {8}\int_0^1 \frac{\log x\log^2\left(\frac {1+x}{1-x}\right)}{x}dx -\frac{\pi^4}{64}. \end{align*}
This gives \begin{align*} J =& -K -\frac{\pi^4}{256}\\ =&-\frac 1 {8}\int_0^1 \frac{\log x\log^2\left(\frac {1+x}{1-x}\right)}{x}dx +\frac{3\pi^4}{256} \end{align*} and \begin{align*} I_2 =& -J +\frac{\pi^4}{64}\\ =&\boxed{\frac 1 8 \int_0^1 \frac{\log x \log^2\left(\frac {1+x}{1-x}\right)}{x}dx +\frac{\pi^4}{256}} \end{align*} as wanted. Thus it holds $$ I = \frac {5\pi^4}{2304} $$ by the previous argument.