I am struggling to solve the following optimization problem:
\begin{align} \max_{ x,y \in [b,1],\frac{p}{x} + \frac{1-p}{y} =a, p \in [0,1]} \left( p x^2 +(1-p) y^2 \right) \end{align} where we also assume that $a>1,b>0$ and $a\le \frac{1}{b}$ (this is needed for the constraint set to be non-empty). I tried solving solving for $p$ in terms of $(x,y)$, things are not exactly simple. If you solve for $p$ from $\frac{p}{x} + \frac{1-p}{y} =a$ one still needs to keep in mind that $p \in [0,1]$.
The maximum is $$ 1+b+b^2-ab(1+b)=1-(a-1)\big(b+b^2\big) $$ and if $\ ab<1\ $ it is achieved at $\ x=b, y=1, \,p=\frac{b(a-1)}{1-b}\ $ and at $\ x=1,y=b,\,p=\frac{1-ba}{1-b}\ $. When $\ ab<1\ $ these are the only values of $\ x,y\ $ and $\ p\ $ that will achieve the maximum. In this case, note that since $\ a>1\ $ and $\ 0<b<\frac{1}{a}\ $, then \begin{align} 1-(a-1)\big(b+b^2\big)&>1-(a-1)\left(\frac{1}{a}+\frac{1}{a^2}\right)\\ &=\frac{1}{a^2}\ , \end{align} so the stationary point identified in Ninad Munshi's comment is not actually a maximum. When $\ ab=1\ $ the only feasible values of $\ x\ $ and $\ y\ $ are $\ x=y=b=\frac{1}{a}\ $, and the objective reduces to $\ \frac{1}{a^2}\ $ for all values of $\ p\in[0,1]\ .$
You can solve the problem by using the equation $\ \frac{p}{x}+$$\,\frac{1-p}{y}=$$\,a\ $ to express $\ p\ $ in terms of $\ y\ $ and $\ x\ $. From this equation you get $\ p=\frac{x(ya-1)}{y-x} \ $, when $\ y\ne x\ $ or $\ x=y=\frac{1}{a}\ $ if $\ y=x\ .$ When $\ y>x\ $ then, since $\ x\ge b>0\ $, to get $\ 0\le p\le1\ $ you will have to have $\ y\ge\frac{1}{a}\ $ and $$ x(ya-1)\le y-x\ , $$ or, equivalently, $$ x\le\frac{1}{a}\ . $$ Conversely, by essentially the same argument, when $\ y<x\ $ you will require $\ y\le\frac{1}{a}\ $ and $\ x\ge\frac{1}{a}\ $.
In fact, without loss of generality, you can restrict $\ x\ $ to the interval $\ \left[b,\frac{1}{a}\right]\ $ and $\ y\ $ to the interval $\ \left[\frac{1}{a},1\right]\ $, because the problem is symmetric with respect to the interchange of $\ x,p\ $ and $\ y,1-p\ $—that is, if $\ x,y,p\ $ is any feasible triple, then $\ x',y',p'\ $ with $\ x'=y,y'=x\ $ and $\ p'=1-p\ $ is also a feasible triple with the same objective value. Thus, substituting $\ p=\frac{x(ya-1)}{y-x}\ $ into the objective function reduces the problem to the following: \begin{align} &\max_{x\in\left[b,\frac{1}{a}\right]}\max_{y\in\left[\frac{1}{a},1\right]}\frac{x^3(ya-1)+y^3(1-ax)}{y-x}\\ =&\max_{x\in\left[b,\frac{1}{a}\right]}\max_{y\in\left[\frac{1}{a},1\right]}y^2+x^2+xy-xya(y+x) \end{align} with the understanding that for $\ y=x=\frac{1}{a}\ $ the value of the objective is to be taken as $\ \frac{1}{a^2}\ $.
Now \begin{align} y^2+x^2+xy-x&ya(y+x)\\ &=(1-ax)y^2+x(1-ax)y+x^2\ , \end{align} and since $\ 1-ax>0\ $ when $\ x<\frac{1}{a}\ $, then $$ \max_{y\in\left[\frac{1}{a},1\right]}(1-ax)y^2+x(1-ax)y+x^2 $$ occurs at $\ y=1\ $ in that case. Therefore \begin{align} \max_{y\in\left[\frac{1}{a},1\right]}&(1-ax)y^2+x(1-ax)y+x^2\\ &=(1-a)x^2+(1-a)x+1\ . \end{align} (It's easy to check that this is still the case when $\ x=\frac{1}{a}\ $.)
Since $\ a>1\ $, the maximum of $\ (1-a)x^2+$$\,(1-a)x+$$\,1\ $ over the interval $\ \left[b,\frac{1}{a}\right]\ $ occurs at the nearest end-point of that interval to the global maximum at $\ x=-\frac{1}{2}\ $, namely at $\ x=b\ $. Thus, the maximum occurs at $\ x=b, y=1\ $, as stated above, and from the above-mentioned symmetry, it also occurs at $\ x=1, y=b\ $.