I'm interested in solving a boundary control problem for an axisymmetric diffusion problem where diffusive fluxes only appear radially. The corresponding problem for a uni-dimensional slab can be summarized as follows:

Let $y_d$ be a desired temperature distribution at the end time $T$ and $y_0$ be the initial temperature of the body. To find a control $u: \Sigma_T \rightarrow \mathbb{R}$ that minimizes the distance of the actual temperature $y(T,\cdot)$ at the end time and the desired temperature $y_d$, we consider the following optimization problem.

$$ \text{min} \,J(y,u) := \frac{1}{2}\int_\Omega(y(T,x)-y_d(x))^2dx + \frac{\alpha}{2}\int_0^T\int_{\partial \Omega} u(t,x)^2 dS(x) dt \\ \text{subject to}\begin{align}\quad y_t-\Delta y&=0\quad\text{on}\;\Omega_T,\\ \frac{\partial y}{\partial \nu}&=\frac{\beta}{\kappa}(u-y) \quad\text{on}\; \Sigma_T \\ y(0,x) &= y_0(x)\quad\text{on}\;\Omega\\ a &\leq u \leq b\quad\text{on}\;\Sigma_T \end{align} $$

To adapt this problem to cylindrical coordinates, does it suffice to specialize the Laplacian (or a more general non-linear diffusion operator) and the cumulative temperature difference (the first term in the objective function $J(y,u)$) as shown below?

$$ \Delta := \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial}{\partial r}\right) $$

$$ \frac{1}{2} \int_0^{2\pi} d\theta\int_0^R (y(T,r) - y_d(r))^2 \,r \,dr = \pi \int_0^R (y(T,r) - y_d(r))^2 \,r \,dr $$

Is there a numerical code or algorithm you would recommend to solve such problems (or a more general variant with non-linear diffusivity, nonlinear Robin boundary conditions, and additional coupled ODEs) in practice? How do I pick the control regularization parameter $\alpha$?

My current plan was to adopt a first discretize, then optimize approach, using the COBYLA algorithm to solve a constrained optimization problem. This seems to be easier than going down the variational path given by optimal control theory of PDEs.

  • $\begingroup$ Are the controls $u$ also independent of $\theta$? If so, then you're really looking at a 1D optimization problem for the (single) value that $u$ takes. If the controls do depend on $\theta$, then your solution $y$ should as well, which would mean that your final expression for the objective is incorrect. $\endgroup$ Aug 12, 2022 at 15:49
  • $\begingroup$ The controls are independent of $\theta$. Imagine the drying process for a product such as dry spaghetti. $y$ is the moisture concentration, and $u$ is the moisture concentration of the air in the dryer. $\endgroup$
    – IPribec
    Aug 13, 2022 at 5:46

1 Answer 1


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 ;)

  • $\begingroup$ The actual problem has non-constant diffusivity and the boundary control is indirect, i.e. the Robin condition uses $y'(u)$ instead of $u$, where $y'$ is a nonlinear function - the moisture sorption isotherm. For the linear case, solutions of the PDE are given by Crank in The Mathematics of Diffusion. I didn't want to complicate my question further with such aspects. The discretization of my controls is "constrained" by what can be built. A series of drying chambers with fixed $u$ values is akin to a C0 finite element discretization in time. Thank you for giving a broad answer, it helps me! $\endgroup$
    – IPribec
    Aug 15, 2022 at 10:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.