Too long to comment:
You might want to try out the following sub-optimal approaches, as the original problem is likely to be in NP-class (citations needed!). I am also going to assume that $A$ has only real entries, for simplicity.
1) Take a look at the paper "Orthogonalization of vectors with minimal adjustment" in Biometrika. The first optimality criterion might be a good surrogate for your cost function. Intuitively, it finds the an orthogonal basis which has least deviation (in terms of dot product) from the a given set of vectors (here columns of $A$).
2) You might also want to investigate the sub-optimality of $Q$ where $A=QR$ (QR decomposition). Changing the order in which Gram-Schmidt is done will yield different QR decompositions, and so you might also want to use this degree of freedom judiciously.
3) Lastly, the greedy way. Intuitively, suppose one were to find a unit vector $q_1$ such that $||q_1A||_1$ is maximized. Post this, one can find a unit vector maximizing $||q_2A||_1$, subject to the additional constraint $q_1q^\top_2=0$. This can be done till one obtains a orthonormal matrix. What remains now is:
$$
\max_{q\in R^{n\times 1}} ||qA||_1~\mbox{subject to}~qq^\top = 1 ~\&~ qC = 0.
$$
Intuitively, a unit vector has highest 1-norm if all its entries are equal in magnitude, or its alignment along on of the vectors $[\pm 1,\cdots, \pm 1]$ is maximum. This can be used in the following iterative routine. Choose any unit vector $q^{(0)}$ such that $q^{(0)}C=0$. Run the following convex routine iteratively till convergence (there is no guarantee of convergence):
$$
q^{(k+1)} = \arg \max_{qA\in R^{n\times 1}} q\left({\mbox{sign}(q^{(k)})}^\top\right) ~\mbox{subject to}~ qq^\top\leq 1 ~\&~ qC=0.
$$
(4) One could also try the following simpler idea. Note that
$\begin{bmatrix} I & X\\X^\top &I \end{bmatrix} \succeq 0 $ is a convex relaxation to orthonormal constraint (using Schur complement). And a the maximum of a linear cost functional with this relaxed constraint will always yeild an optima at its boundary, which is the set of orthonormal matrix. Using the same logic as in point (3), one can write the code below:
import numpy as np
import cvxpy as cvx
X = cvx.Variable((N,N))
Q,R = np.linalg.qr(np.random.randn(N,N))
Q = np.sign(Q@A)
count = 0
while(count<=20):
constraints = [cvx.vstack((cvx.hstack((np.eye(N),X)),cvx.hstack((X.T,np.eye(N))))) >> 0]
prob = cvx.Problem(cvx.Maximize(cvx.trace((X @ A)@(Q.T))), constraints)
prob.solve()
Q = np.sign((X.value)@A)
count = count + 1
print(np.sum(np.abs(X.value @ A)))
Z,R = np.linalg.qr(np.random.randn(N,N))
print(np.sum(np.abs(Z@A)))
Hope this helps.