I can reduce this discrete optimization problem to a continuous optimization problem in many cases by using my own techniques, but the continuous optimization problem may return a higher value than the discrete optimization problem in some cases. Sometimes the continuous optimization problem finds a solution that does not exist in the discrete case.

The OP has observed that if R,S are matrices, then $\text{Tr}(RS)=\sum(R\odot S^T)$ and therefore $\text{Tr}(ADBD)=\sum(D(A\odot B^T)D)$.

The problem of maximizing $\sum(DCD)$ where $D$ is a 0-1 diagonal matrix with $\text{Tr}(D)=m$ is NP-complete since the clique problem in graph theory is a special case of this problem. But as with many NP-complete problems, one can develop heuristic algorithms for maximizing $\sum(DCD)$ such as hill climbing algorithms or evolutionary algorithms.

The expression $\max_D\sum(DCD)$ is a limit case of the expression $\max_D\rho(DCD)$ where $\rho$ denotes the spectral radius. If $N_n$ is the $n\times n$-matrix where each entry in $N_n$ is $1$, and $C$ is a real matrix, then $\sum(C)=\lim_{t\rightarrow 0}\frac{\rho(nt\cdot C-N_n)-n}{t}=\frac{\partial}{\partial t}\rho(nT\cdot C-N_n)|_{t=0}$ (this limit can be derived from the formula for the gradient of the eigenvalues of a matrix), and we have a similar expression for $\max_{D}\sum(DCD)$ in terms of the spectral radius.

Suppose that $K$ is either the field of real numbers or complex numbers (we will need to make modifications to generalize everything to quaternionic matrices) and $(A_1,\dots,A_r),(B_1,\dots,B_r)\in M_n(K)^r$. Then define the $L_2$-spectral radius of $(A_1,\dots,A_r)$ by $\rho_2(A_1,\dots,A_r)=\rho(A_1\oplus\overline{A_1}+\dots+A_r\oplus\overline{A_r})$. If $(A_1,\dots,A_r)\in M_m(K),(B_1,\dots,B_r)\in M_n(K)$, then define the $L_2$-spectral radius similarity between $(A_1,\dots,A_r)$ and $(B_1,\dots,B_r)$ is the quantity $$\|(A_1,\dots,A_r)\simeq(B_1,\dots,B_r)\|_2=\frac{\rho(A_1\oplus\overline{B_1}+\dots+A_r\oplus\overline{B_r})}{\rho_2(A_1,\dots,A_r)\rho_2(B_1,\dots,B_r)}.$$

Define the $L_{2,d}$-spectral radius of $A_1,\dots,A_r$ by setting $$\rho_{2,d}^K(A_1,\dots,A_r)=\max\frac{\rho(A_1\oplus\overline{X_1}+\dots+A_r\oplus \overline{X_r})}{\rho_2(X_1,\dots,X_r)}:X_1,\dots,X_r\in M_d(K)\}.$$

If $(X_1,\dots,X_r)\in M_d(K)^n$ locally maximizes the quantity $$\frac{\rho(A_1\oplus \overline{X_1}+\dots+A_r\oplus\overline{X_r})}{\rho_2(X_1,\dots,X_r)},$$ then we say that $(X_1,\dots,X_r)$ is an $L_2$-spectral radius dimensionality reduction (LSRDR) of $(A_1,\dots,A_r)$. We observe that $(X_1,\dots,X_r)$ is an LSRDR of $(A_1,\dots,A_r)$ precisely when $(X_1,\dots,X_r)$ maximizes the similarity $\|(A_1,\dots,A_r)\simeq(X_1,\dots,X_r)\|$.

Empirical observation: Suppose that $(X_1,\dots,X_r)\in M_d(K)^r$ is an LSRDR of $(A_1,\dots,A_r)$. Then we can typically find matrices $R,S$ and some constant $\lambda$ where $X_j=\lambda RA_jS$ for all $j$ and where $RS$ is the identity matrix, and $SR$ is a (typically non-orthogonal) projection. In this case, let $P=SR$. The projection matrix $P$ is often unique in the sense that if we compute the LSRDR again and compute another projection matrix $P_1$, then we would get $P=P_1.$ Observe that $\rho_2(X_1,\dots,X_r)=\rho_2(PA_1P,\dots,PA_rP)$ and $\rho(A_1\otimes\overline{X_1}+\dots+A_r\otimes\overline{X_r})=\rho(A_1\otimes\overline{PA_1P}+\dots+A_r\otimes\overline{PA_rP})$.

Given $A_1,\dots,A_r\in M_n(K),X_1,\dots,X_r\in M_d(K)$, define a superoperator $$\Gamma(A_1,\dots,A_r;X_1,\dots,X_r):M_{n,d}(K)\rightarrow M_{n,d}(K)$$ by letting $$\Gamma(A_1,\dots,A_r;X_1,\dots,X_r)(X)=\sum_{k=1}^rA_kXX_r^*,$$ and set $$\Phi(A_1,\dots,A_r)=\Gamma(A_1,\dots,A_r;A_1,\dots,A_r).$$

Suppose now that $C=(c_{i,j})_{i,j}$ is a matrix with non-negative real entries. Let $a=(a_{i,j})_{i,j}$ be a matrix with $|a_{i,j}|^2=c_{i,j}$ for all $i,j$. For each $i,j$, let $A_{i,j}$ be the matrix where each entry except for the $i,j$ entry is $0$, but where the $(i,j)$ entry of $A_{i,j}$ is $a_{i,j}$. Suppose now that $(X_{i,j})\in M_d(K)^{n\times n}$ is an LSRDR of $(A_{i,j})_{i,j}$. Then the projection matrix $P$ is often a diagonal projection matrix where the diagonal entries are either $0$ or $1$.

If $v$ is a vector, then let $\text{Diag}(v)$ denote the square matrix whose diagonal entries are $v$, and if $V$ is a matrix, then let $\text{Diag}(V)$ denote the vector whose entries are the diagonal entries in $V$.

$\Gamma((A_{i,j})_{i,j},(PA_{i,j}P)_{i,j})=\Phi((PA_{i,j}P)_{i,j})$, and $\Phi((PA_{i,j}P)_{i,j})(X)=\text{Diag}(C\cdot\text{Diag}(X))$ for all matrices $X$. This means that $$\rho(\Gamma((A_{i,j})_{i,j},(PA_{i,j}P)_{i,j}))=\rho(\Phi((PA_{i,j}P)_{i,j})=\rho(PCP),$$ so

$$\max\{\rho\big(\text{Diag}(D)\cdot C\cdot\text{Diag}(D)\big):D\in\{0,1\}^n,\sum D=d\}:=\rho_{2,d}((A_{i,j})_{i,j})^{1/2}$$ if everything works well.

I may add details to this answer about how we can write $$\max\{\rho\big(\text{Diag}(D)\cdot C\cdot\text{Diag}(D)\big):D\in\{0,1\}^n,\sum D=d\}$$ as a continuous optimization problem without needing to use the spectral radius.

You can find more information about LSRDRs at my cryptocurrency site circcashcore.org/blog (I use LSRDRs and similar machine learning algorithms to evaluate the cryptographic security of block ciphers).