MLE for a modified German tank problem Suppose I have a distribution $U(0,a)$ where $a$ is unknown and we are interested to estimate it. Someone who has access to $n$ samples $\mu_i$ of this distribution instead decides to create a Gaussian distribution with mean $\mu_i$ and standard deviation $\sigma_i$ and draws 1 sample from each of $x_i\sim N(\mu_i, \sigma_i)$ and provides us this $x_i, \sigma_i$. Based on this we would like to estimate distribution of $a$.
Since this problem is quite similar to estimating the distribution of maximum, my approach was to write the CDF of $a$ as the product of individual of CDFs : $$ \phi(a) = \prod_{i=1}^n \phi(a, x_i, \sigma_i)$$, but I know this doesn’t capture the fact that there is a cutoff at $a$. Any ideas on how to improve on this?

Given your specified distributions, the marginal density of $X$ is:
$$\begin{align}
f_X(x) 
&= \int \limits_0^a \text{N}(x |\mu,\sigma^2) \cdot \text{U}(\mu|0,a) \ d \mu \\[6pt]
&= \frac{1}{a} \int \limits_0^a \frac{1}{\sigma \sqrt{2 \pi}} 
 \exp \bigg( -\frac{1}{2} \bigg( \frac{x-\mu}{\sigma} \bigg)^2 \bigg) \ d \mu \\[6pt]
&= -\frac{1}{a} \int \limits_{x/\sigma}^{(x-a)/\sigma} \frac{1}{\sqrt{2 \pi}} 
 \exp \bigg( -\frac{1}{2} \bigg( \frac{x-\mu}{\sigma} \bigg)^2 \bigg) \ d ((x-\mu)/\sigma) \\[6pt]
&= \frac{1}{a} \bigg[ \Phi \Big( \frac{x}{\sigma} \Big) - \Phi \Big( \frac{x-a}{\sigma} \Big) \bigg]. \\[6pt]
\end{align}$$
So your log-likelihood function should be:
$$\ell_{\mathbf{x}}(a)
= - n \log(a) + \sum_{i=1}^n \log \bigg[ \Phi \Big( \frac{x_i}{\sigma_i} \Big) - \Phi \Big( \frac{x_i-a}{\sigma_i} \Big) \bigg].$$
Note that there is no upper bound on $a$ in this problem because the $x_i$ values all have support over the whole real line irrespective of the value of $a$.  To find the MLE, your score function and information function are:
$$\begin{align}
s_{\mathbf{x}}(a)
\equiv \frac{d \ell_{\mathbf{x}}}{da} (a)
&= - \frac{n}{a} + \sum_{i=1}^n \frac{\phi((x_i-a)/\sigma_i)}{\sigma_i H_\mathbf{x}(a)}, \\[12pt]
I_{\mathbf{x}}(a)
\equiv -\frac{d^2 \ell_{\mathbf{x}}}{da^2} (a)
&= - \frac{n}{a^2} - \sum_{i=1}^n (x_i-a) \cdot \frac{\phi((x_i-a)/\sigma_i)}{\sigma_i^3 H_\mathbf{x}(a)} + \sum_{i=1}^n \frac{\phi((x_i-a)/\sigma_i)^2}{\sigma_i^2 H_\mathbf{x}(a)^2}, \\[12pt]
\end{align}$$
where $H_\mathbf{x}(a) \equiv \Phi(x_i/\sigma_i) - \Phi((x_i-a)/\sigma_i)$.  Taking $s_{\mathbf{x}}(\hat{a}) = 0$ and solving numerically for $\hat{a}$ will yield the MLE for $a$.  You should bear in mind that the MLE tends to be a biased estimator for problems involving estimating a maximum-bound parameter, but it is a common form of estimation in general.
Using the laws of iterated expectation and variance, it is simple to show that $\mathbb{E}(X) = \tfrac{1}{2}a$ and $\mathbb{V}(X) = \tfrac{1}{12} a^2 + \sigma^2$.  Consequently, under some relatively weak conditions on the values $\sigma_i$ (ensuring that we don't have a finite subset of values that "dominate" in the limit) we have $2\bar{x}_n \rightarrow a$ as $n \rightarrow \infty$.  Consequently, for large $n$ we ought to have $\hat{a} \approx 2\bar{x}_n$.  (But be careful with this approximation; in some cases we may have $\bar{x}_n \leqslant 0$.)

Implementation in R: We can implement the numeric solution to the MLE for this problem using the following R function.  This function takes in the vectors $\mathbf{x}$ and $\boldsymbol{\sigma}$ and returns the MLE $\hat{a}$ and the mean-maximum-log-likelihood $\ell_{\mathbf{x}}(\hat{a})/n$.  (Optimisation is done here using the nlm function rather than using a root-finding algorithm on the score function.)  The function also has an option to return the convergence code from the optimisation to check that the numerical maximisation worked properly.
MLE.modified.tank <- function(x, sigma = rep(1, length(x)), nlm.code = FALSE) {
  
  #Check inputs
  if (!is.numeric(x))             stop('Error: Input x should be numeric')
  if (!is.numeric(sigma))         stop('Error: Input sigma should be numeric')
  if (length(x) != length(sigma)) stop('Error: Inputs x and sigma should be the same length')
  if (min(sigma) <= 0)            stop('Error: Input sigma should only have positive values')   
  
  #Set objective function using parameterisation p = log(a)
  NEGLOGLIKE <- function(p) {
    n <- length(x)
    LOG1 <- pnorm(x/sigma, log = TRUE)
    LOG2 <- pnorm((x-exp(p))/sigma, log = TRUE)
    LOGDIFF <- LOG1 + VGAM::log1mexp(LOG1 - LOG2)
    n*p - sum(LOGDIFF) }
  
  #Find the MLE
  P0    <- log(max(2*mean(x), 1e-6))
  NLM   <- nlm(f = NEGLOGLIKE, p = P0)
  MLE.a <- exp(NLM$estimate)
  MMAX  <- -NLM$minimum/n
  CODE  <- NLM$code
  
  #Create the output
  OUT <- data.frame(MLE.a = MLE.a, mean.max.loglike = MMAX)
  if (nlm.code) { OUT$nlm.code <- CODE }
  rownames(OUT) <- ''
  
  #Return the output
  OUT }

We can implement this function to see if it estimates well.  In the code below we take a set of mock data generated from your model (with known parameter) and we estimate via MLE.  With $n=200$ data points this gives a good estimate of the true parameter.
#Set parameters
a <- 20
n <- 200
sigma <- rgamma(n, shape = 6, scale = 0.4)

#Generate some mock data
set.seed(1)
m <- runif(n, min = 0, max = a)
x <- m + sigma*rnorm(n)

#Compute the MLE of a
MLE <- MLE.modified.tank(x, sigma, nlm.code = TRUE)

#Plot the data and MLE
plot(sigma, x, 
     main = '(Blue line shows true parameter - Red line shows MLE)',
     xlab = 'SDs of observed values',
     ylab = 'Observed values')
abline(h = a, lty = 2, col = 'blue')
abline(h = MLE$MLE.a, col = 'red')

#Display the MLE
MLE 

   MLE.a mean.max.loglike nlm.code
 19.7711        -3.178023        1


