I did the case $p=1$ here. The generalization to higher $p$ may involve higher-order derivatives as follows:

$$\begin{align}K_p &= \int_0^{\pi/2} dx \frac{x^{2 p}}{1+\cos^2{x}} = \frac1{2^{4 p-1}} \int_{-\pi}^{\pi} dy \frac{y^{2 p}}{3+\cos{y}} \end{align}$$

So define, as before,

$$J(a) = \int_{-\pi}^{\pi} dy \frac{e^{i a y}}{3+\cos{y}} $$

Then

$$K_p = \frac{(-1)^p}{2^{4 p-1}} \left [\frac{d^{2 p}}{da^{2 p}} J(a) \right ]_{a=0}$$

In the above reference, I derive an expression for $J(a)$:

$$J(a) = \frac{\pi}{\sqrt{2}} \cos{\pi a} \: (3-2 \sqrt{2})^a – 4 \sin{\pi a}\: PV \int_0^1 dx \frac{x^a}{x^2-6 x+1}$$

where $PV$ denotes the Cauchy principal value. Thus, the above outlines a procedure for expressing the integral $K_p$ in terms of powers of $\log{(3-2 \sqrt{2})}$ and $\operatorname{Li}_{2 k}(3-2 \sqrt{2})$ for $k=1$ to $p$. I hope I can derive a general expression, but for the time being, this outline should suffice for arbitrary positive integer values of $p$.

**Evaluating the derivatives**

I have found a formalism in which I can derive a more explicit expression for the $K_p$. Consider the following function:

$$f(a) = \frac{\pi}{\sqrt{2}} \cos{\pi a} \: (3-2 \sqrt{2})^a – 4 \sin{\pi a} \, x^a$$

We seek to take $p$ successive second derivatives of $f$ with respect to $a$. Note how both terms of $f$ are a trig function times a number $q$ raised to the $a$th power. It turns out that

$$\frac{d^2}{da^2} \left ( \begin{array} & \sin{\pi a} \, q^a \\ \cos{\pi a} \, q^a \end{array}\right ) = \left ( \begin{array} & \log^2{q}-\pi^2 & 2 \pi \log{q} \\ -2 \pi \log{q} & \log^2{q}-\pi^2 \end{array} \right ) \left ( \begin{array} & \sin{\pi a} \, q^a \\ \cos{\pi a} \, q^a \end{array}\right )$$

We can express this in terms of a rotation matrix:

$$\frac{d^2}{da^2} \left ( \begin{array} & \sin{\pi a} \, q^a \\ \cos{\pi a} \, q^a \end{array}\right ) = \left ( \log^2{q}+\pi^2\right )\left ( \begin{array} & \cos{\theta} & \sin{\theta} \\ -\sin{\theta} & \cos{\theta} \end{array} \right ) \left ( \begin{array} & \sin{\pi a} \, q^a \\ \cos{\pi a} \, q^a \end{array}\right )$$

where $\cos{\theta} = (\log^2{q}-\pi^2)/(\log^2{q}+\pi^2)$ and $\sin{\theta} = 2 \pi \log{q}/(\log^2{q}+\pi^2)$. We then have an explicit expression for the $2 p$th derivative of the terms of $f$:

$$\begin{align}\left [\frac{d^{2 p}}{da^{2 p}} \left ( \begin{array} & \sin{\pi a} \, q^a \\ \cos{\pi a} \, q^a \end{array}\right )\right ]_{a=0} &= \left ( \log^2{q}+\pi^2\right )^p\left ( \begin{array} & \cos{p \theta} & \sin{p \theta} \\ -\sin{p \theta} & \cos{p \theta} \end{array} \right ) \left ( \begin{array} & 0 \\ 1 \end{array}\right )\\ &= \left ( \log^2{q}+\pi^2\right )^p \left ( \begin{array} & \sin{p \theta} \\ \cos{p \theta} \end{array}\right ) \end{align}$$

The above elements are simply Chebyshev polynomials of the second and first kind, respectively, in $\cos{\theta}$. As an example, let’s do the first few values of $p$:

$p=1$:

$$\left [\frac{d^{2}}{da^{2}} f(a) \right ]_{a=0} = \frac{\pi}{\sqrt{2}} \left (\log^2{(3-2 \sqrt{2})}-\pi^2 \right ) – 8 \pi \log{x}$$

which agrees with the result in the linked solution.

$p=2$:

$$\left [\frac{d^{4}}{da^{4}} f(a) \right ]_{a=0} = \frac{\pi}{\sqrt{2}} \left [ \left (\log^2{(3-2 \sqrt{2})}-\pi^2 \right )^2-4 \pi^2 \log^2{(3-2 \sqrt{2})}\right ] – 16 \pi \log{x} \left (\log^2{x}-\pi^2 \right )$$

$p=3$:

$$\left [\frac{d^{6}}{da^{6}} f(a) \right ]_{a=0} = \frac{\pi}{\sqrt{2}} \left [\left (\log^2{(3-2 \sqrt{2})}-\pi^2 \right )^3-12 \pi^2 \log^2{(3-2 \sqrt{2})} \left (\log^2{(3-2 \sqrt{2})}-\pi^2 \right ) \right ]

\\ – 4 \left [6 \pi (\log^2{x}-\pi^2)^2 \log{x} - 8 \pi^3 \log^3{x} \right ]$$

and so on. Once this is done, the $\log{x}$ terms need to thrown into the integral above, i.e., we need to compute

$$H_k = PV \int_0^1 dx \frac{\log^{2 k-1}{x}}{x^2-6 x+1} $$

for $k = 1,\ldots,p$. Each $H_k$ will be expressible in terms of a $\operatorname{Li}_{2 k}(3-2 \sqrt{2})$ plus some other terms. If I have more time, I will post an explicit expression for $H_k$ and finally get an explicit expression for the original integral. But for now, this represents a lot of progress.

**Integration**

OK, it is in fact practical and relatively straightforward, albeit a huge mess, to evaluate the $H_k$. The basis for the evaluation lies in the expression of $H_k$ as

$$H_k = \frac1{4 \sqrt{2}} \left [\int_0^1 dx \frac{\log^{2 k-1}{x}}{x-(3+2 \sqrt{2})} - PV \int_0^1 dx \frac{\log^{2 k-1}{x}}{x-(3-2 \sqrt{2})} \right ]$$

The first integral is straightforward:

$$\int_0^1 dx \frac{\log^{2 k-1}{x}}{x-(3+2 \sqrt{2})} = (2k-1)! \operatorname{Li}_{2 k}(3-2 \sqrt{2})$$

The 2nd one, however, is somewhat challenging. Note that the principal value carries into this integral because there is a bona fide singularity, but recall that the principal value was not just are up to avoid this singularity. Rather, it was a consequence of the way we treated the pole on the branch cut.

Anyway, the second integral is given by

$$PV \int_0^1 dx \frac{\log^{2 k-1}{x}}{x-(3-2 \sqrt{2})} = – i \pi \log^{2 k-1}(3-2 \sqrt{2}) + (2 k-1)! \operatorname{Li}_{2 k}(3 + 2 \sqrt{2})$$

Now, it turns out that the we may actually determine a value for the polylog of a positive real number greater than $1$ by a general inversion relation:

$$\operatorname{Li}_{2 k}(3 – 2 \sqrt{2}) + \operatorname{Li}_{2 k}(3 + 2 \sqrt{2}) = (-1)^{k+1} \frac{(2 \pi)^{2 k}}{(2 k)!} B_{2 k} \left (1-i\frac1{2 \pi} \log{(3+2 \sqrt{2})} \right ) $$

where $B_{2 k}$ is the $(2 k)$th Bernoulli polynomial. The leads to a formulation for evaluating each $H_k$. Note that the imaginary pieces should cancel out. For example,

$$4 \sqrt{2} H_1 = 2 \operatorname{Li}_2(3-2 \sqrt{2}) – \frac{\pi^2}{3} + \frac12 \log^2{(3-2\sqrt{2})}$$

$$4 \sqrt{2} H_2 = 2 \operatorname{Li}_4(3-2 \sqrt{2}) – \frac{2 \pi^4}{15} – \pi^2 \log^2{(3-2\sqrt{2})} +\frac14 \log^4{(3-2 \sqrt{2})}$$

$$4 \sqrt{2} H_3 = 2 \operatorname{Li}_6(3-2 \sqrt{2}) – \frac{16 \pi^6}{63} – \frac{4 \pi^4}{3} \log^2{(3-2\sqrt{2})} – \frac{5 \pi^2}{3} \log^4{(3-2 \sqrt{2})}+\frac{1}{6} \log^6{(3-2 \sqrt{2})}$$