Consider the problem: "one dice has three faces $1,2,3$ with probabilities $p_1,p_2,p_3$ and it is thrown many times with a mean value $a$. What are the values for $p_1,p_2,p_3$ which maximizes the uncertainty consistent with the mean value $a$?"
Here I believe uncertainty should be understood as standard deviation.
So, the experiment for a single such dice being thrown has outcomes described by the random variable $X$ taking values on $\{1,2,3\}$ with probabilites $p_1,p_2,p_3$ respectively.
I believe that it all boils down to maximizing
$$\sigma^2(p_1,p_2,p_3)=E[X^2](p_1,p_2,p_3)-E[X]^2(p_1,p_2,p_3)$$
subject to the constraints
$$g(p_1,p_2,p_3)=p_1+p_2+p_3=1\quad h(p_1,p_2,p_3)=E[X](p_1,p_2,p_3)=p_1+2p_2+3p_3=a.$$
Thus it seems we should use Lagrange multipliers and solve
$$\nabla \sigma^2=\lambda \nabla g+\mu\nabla h,\\ g(p_1,p_2,p_3)=1, \\ h(p_1,p_2,p_3)=a.$$
Furthermore, it is clear that we have
$$\dfrac{\partial}{\partial p_i}E[f(X)]=\dfrac{\partial}{\partial p_i}\sum_j f(x_j)p_j=f(x_i).$$
so that we have
$$\nabla \sigma^2(p_1,p_2,p_3)=(1-2E[X],4-4E[X],9-6E[X]).$$
Thus we have the system of equations
$$\begin{cases}1-2E[X]&=\lambda+\mu \\ 4-4E[X]&=\lambda+2\mu\\ 9-6E[X]&=\lambda+3\mu\\ p_1+2p_2+3p_3&= a \\ p_1+p_2+p_3&= 1.\end{cases}$$
The fourth equation imposes $E[X]=a$. Hence substituting it on the first and second we can find $\mu = 3-2a$ and $\lambda = -2$.
But now this ought to be wrong. If I substitute this on the third equation I get $$9-6E[X]=-2+9-6a\Longleftrightarrow 6E[X]=6a+2$$
which is incompatible with the fourth equation which tells that $E[X]=a$.
So this method doesn't work. What am I doing wrong here? What is the right way to solve this problem?
In this particular example of three variables $p_i$, $i=1,2,3$, we can solve this explicitly as a linear problem. More general cases where $i=1,2,\ldots,n$ most likely will not succumb to such a simple approach.
We assume the scenario as given above, in which a random variable $X$ takes values $1,\,2,\,3$ with associated probabilities $p_1,\,p_2,\,p_3$.
In your question you wish to consider the probabilities such that:
We note that since variances of independent samples are linear, we have
$$ \text{Var}(X_1 + \cdots + X_N) = N\text{Var}(X_1),$$
from which it follows that it suffices to answer the question: find the probabilities such that
That is: we can solve the question by considering just a single $X$; henceforth we drop the dependence on $N$, since this is an arbitrary constant, so just ask to find the $X$ with mean $\alpha$ and maximal variance.
We want that the mean is equal to some given value $\alpha \in [1,3]$ (noting that it is not possible for $\alpha$ to be outside of this interval, given that $X \in \{1,2,3\}$). Then we have
$$\mathbf E[X] = p_1 + 2p_2 + 3p_3 = \alpha,$$
and moreover we have the following constraints
\begin{align*} p_1+ p_2 +p_3 & = 1, \\ 0 \leq p_1 & \leq 1 ,\\ 0 \leq p_2 & \leq 1 ,\\ 0 \leq p_3 & \leq 1. \\ \end{align*}
Using the first constraint, and the expression for the expectation above we see that both $p_2$ and $p_3$ can be expressed in terms of $p_1,\, \alpha$
\begin{align*} p_2 & = (3-\alpha) - 2p_1, \\ p_3 & = (\alpha -2) + p_1. \end{align*}
Not all combinations $(p_1,\,p_2,\,p_3)$ generated in this way form valid probability distributions: for this we will rely on the remaining constraints: $0 \leq p_i \leq 1$, $i=1,2,3$. Reformulating each of these in terms of $p_1$ leads to
\begin{align*} 0 & \leq p_1 \leq 1. \\ 1 - \frac12\alpha & \leq p_1 \leq \frac32 - \frac12\alpha, \\ 2 - \alpha & \leq p_1 \leq 3 - \alpha. \end{align*}
Depending on the value of $\alpha$ many of these constraints will be redundant (eg. for the second inequality, if $\alpha = 3$ then clearly the left hand side is less than $0$, which is already ruled out by the first constraint).
But most importantly, we have found a description of the system using only a single variable (since $\alpha$ is constant), and with the single (non-linear) constraint
$$\max \left( 0, \,1-\frac12 \alpha, \, 2- \alpha \right) \leq p_1 \leq \min\left(1, \, \frac32 - \frac12 \alpha, 3 - \alpha \right).$$
In fact on checking cases, we note that this simplifies to $$\max \left( 0, \, 2- \alpha \right) \leq p_1 \leq \frac32 - \frac12 \alpha.$$
Turning to the uncertainty, which you describe as being the standard deviation, we note that maximizing this is equivalent to maximizing variance; moreover the variance is given to be
\begin{align*} \text{Var}(X) & = \mathbf E[X^2] - \mathbf E[X]^2 \\ & = \mathbf E[X^2] - \alpha^2 \\ & = \left(p_1 + 4 p_2 + 9p_3\right) - \alpha^2 \\ & = p_1 + 4 \big( (3-\alpha) -2p_1\big) + 9 \big( (\alpha -2 ) + p_1 \big) - \alpha^2\\ & = 2p_1 - 6 + 5\alpha - \alpha^2. \end{align*} We see immediately without taking derivatives that the variance is linear (and increasing) in $p_1$, and therefore has its maximum at the largest feasible value of $p_1$.
That is, the maximal variance random variable $X \in \{1,2,3\}$ with mean $\mathbf E[X] = \alpha$ is determined by:
$$p_1 = \frac32 - \frac12 \alpha, \qquad p_2 = 0,\qquad p_3 = \frac12 \alpha - \frac12. $$ In particular we see that for any $\alpha \in [1,3]$ the solution is always has $p_2 = 0$.
In terms of standard distributions, we see that for fixed $\alpha \in[1,3]$ the maximal variance $X$ is given by
$$X = 1 + 2 Y, \qquad Y \sim \text{Ber}\left( \frac12(\alpha-1) \right).$$
As a final addition, though not strictly a part of the question, we note that from this we can also see that the distribution with maximal variance on $\{1,2,3\}$, is given by the case $X = 1 + 2\text{Ber}(1/2)$.
To see this we note that if $Y \sim \text{Ber}(q)$ for some $0 \leq q \leq 1$, then $X = 1 + 2 Y$ has the same variance as $Y$, and moreover this is
$$\text{Var}(X) = \text{Var}(Y) = q(1-q),$$
and in particular this is maximized at $q = 1/2$, which is equivalent to $\alpha =2$.