There are a couple of questions here, some of which pertain specifically to the geometry and input data for your problem and some have more to do with PDE-constrained optimization in general.
Some of my answers for the latter are more my opinion and some are more fact-based, so I'll try to distinguish between them.
Is your objective correct?
Yes, using the model-data misfit functional
$$J(y) = \pi\int_0^R(y(T, r) - y_d(r))^2\,r\,dr$$
will be correct assuming that your data really are just radially-dependent and (this is the important part) that you can apply some numerical method that will be guaranteed to give you radially-symmetric solutions.
Is your expression for the Laplace operator correct?
Your specific problem has constant coefficients and the boundary data $u$ are constant in space.
As long as $u$ is continuous in time then there ought to be a strong solution by the usual regularity arguments.
So the form of the Laplace operator that you've written down is correct for this problem instance, but it could easily be wrong with relatively small changes in the assumptions.
For example, suppose you lose the radial symmetry or have a non-constant thermal conductivity.
You might then have to think about spatially discontinuous values for the controls or the conductivity, and now we're firmly in weak solution territory.
Is there a numerical code or algorithm you would recommend?
The PDE in your problem has an analytical solution expressed as a series in $\{J_0(\alpha_nr)\}$ where $J_0$ is the 0th order Bessel function of the first kind and $\alpha_n$ need to be chosen to satisfy the boundary conditions.
(Admittedly I'd have to dust off my copy of the Habermann PDE book to fill in all the details here.)
This gives you some options that aren't available for general problems.
If you're able to calculate the coefficients of $y_d$ in terms of this Bessel function series then this can greatly simplify calculating that part of the objective.
At least for this problem, I'd try to go as far as you can using the analytical solution for the PDE so that you can focus on what to do with the controls.
How to pick the regularization parameter?
This partly depends on the application.
For state and parameter estimation, this is often dictated by the measurement errors and people use things like the L-curve or the discrepancy principle.
For more engineering or control theory type applications, this is more application dependent.
It's more likely that you know that
$$\int_0^Tu^2\,dt \le TU^2$$
for some $U$ and the parameter $\alpha$ is merely a scalar Lagrange multiplier included in the problem to enforce that inequality constraint.
How will you discretize the controls?
The controls are constant in space, so this is really just a question of how to discretize them in time.
You could use piecewise polynomials of more or less arbitrary degree; you could divide the time interval into however many sub-intervals you like; and you could impose whatever degree of continuity between sub-intervals you like, so you might take $u$ to be discontinuous, continuous, continuously differentiable, etc.
This would be a finite element-type discretization.
Alternatively, you could use spectral methods -- Legendre polynomials or Fourier series.
The choice of discretization is going to depend on how regular you expect $u$ to be.
Your regularization functional is just the square norm of $u$, which doesn't explicitly penalize lack of smoothness.
So you'd probably want to use a finite element discretization with a low or no degree of continuity.
If your regularization functional also penalized the squared derivative or second derivative of $u$, then you might instead think about a Hermite spline or spectral discretization.
The box constraint $a \le u \le b$ also effectively restricts you to piecewise constant or piecewise linear basis functions.
The relationship between the maximum and minimum values of the function and of its coefficients is much more of a pain to evaluate at higher order.
You might want to check this but I'm reasonably sure that you can't expect a rate of convergence faster than $\mathscr{O}(\delta t)$ for discretizations of variational inequalities anyway, so there's no point attempting to go to higher degree.
Can you expect success from a general constrained optimization method?
In my experience, no, you cannot.
The natural vector space in which to look for solutions of PDE or of PDE-constrained optimization problems is a Sobolev space.
These vector spaces are infinite-dimensional and that means weird things start to happen.
The value of the derivative $\partial J/\partial u$ of your objective lives in the dual $Q^*$ of the vector space $Q$ where the controls live.
Now in finite dimensions we grow accustomed to thinking that $Q^*$ is identical to $Q$.
In infinite dimensions this is no longer true, and trying to naively apply common optimization algorithms like gradient descent leads to nonsensical statements like adding an element of $Q^*$ to an element of $Q$.
The convergence rate of the resulting algorithm deteriorates with problem size.
Now this is partly my opinion, but I find the results of using these algorithms without accounting for the function space to be almost completely unreliable.
Yes, the objective does decrease from one iteration to the next, but because of the very slow convergence rate I find it very difficult to tell when you should stop iterating.
Compounding this issue is the fact that, by and large, people who work on optimization methods only ever think about finite-dimensional problems.
The machine learning craze also has a lot of people fixated on first-order methods and those types of algorithms often don't make sense in function spaces, although some do, for example mirror descent.
One remedy is to use second-order methods like Gauss-Newton rather than first-order methods like gradient descent.
Regarding the right map between the dual and primal spaces, ROL is one of the few software packages I know of that provides the right hooks to deal with this issue.
Also, in future, you'll get shorter and better responses if you stick to one question per Question ;)
