Let $T$ denote a $J$-component categorical random variable with pmf $$ \mathsf P(T=t_j)=p_j,\quad j=1,2,\dots,J, $$ where $t_j\in[0,t_\max]$, $t_\max>0$.

I came across a problem that seeks to determine the pmf for $T$, i.e., the optimal probabilities $p_j$ and support $t_j$, that maximizes $$ f(\mathbf p,\mathbf t)=(\mathsf ET)(\mathsf ET+\Phi\mathsf{Var}T), $$ for a specified choice of $\Phi>0$.

I started with a numerical implementation of the optimization problem to see what kinds of observations could be made (the code below minimizes $-f$ as opposed to maximizing $f$). What's interesting, is that for certain choices of $t_\max$ and $\Phi$ the maximum value of $f$ appears to be independent of $J$. Consider the code below using $t_\max=5$ and $\Phi=5$. For $J=2,3,4,5,6$ (possibly for all $J\geq 2$), the maximum value for $f$ is constant and the optimal $t_j$ are all equal to $0$ or $t_\max$. If any of the $t_j$ are not equal to these boundary values, they are assigned a probability of zero. Can anyone present a reason for this behavior? Note that $\mathsf ET$ is maximized when $t_j=t_\max$ for all $j$, while $\mathsf{Var}T$ is maximized when all the $t_j$ are at the boundaries $t_j=0$ and $t_j=t_\max$, so this seems like a clue.

% optimization parameters
    global J
    J = 2;     % number of samples
    Phi = 5;   % flux (e-/s)
    t_max = 5; % maximum allowable exposure time (s)

% minimize generalized var. w.r.t. sample weights and exp. times
    t_bar = @(w,t) sum(w.*t);
    t_hat = @(w,t) sum(w.*(t-t_bar(w,t)).^2);
    obj = @(x) -t_bar(x(1:J),x((J+1):(2*J)))...
        *(t_bar(x(1:J),x((J+1):(2*J)))+Phi*t_hat(x(1:J),x((J+1):(2*J))));

    w0 = rand(1,J);
    w0 = sort(w0)/sum(w0);
    t0 = sort(rand(1,J)*t_max);
    x0 = [w0 t0];
    lb = [zeros(1,J) zeros(1,J)];
    ub = [ones(1,J) t_max*ones(1,J)];

    [sol,val] = fmincon(obj,x0,[],[],[],[],lb,ub,@norm);

% value of obj. function at minimum
    val
    
% optimal sample size weights and exp. times
    w_opt = sol(1:J)
    t_opt = sol((J+1):(2*J))

% constraint function for optimization above
    function [c,ceq] = norm(x)
    global J
    c = [];
    ceq = sum(x(1:J))-1;
    end