The problem is given by:
$$ \arg \min_{X \in \mathcal{T} } \frac{1}{2} {\left\| X B - A \right\|}_{F}^{2} $$
Where $ \mathcal{T} $ is the set of Lower Triangular Matrices.
The set $ \mathcal{T} $ is a Convex Set.
Moreover, the orthogonal projection onto the set of a given matrix $ Y \in \mathbb{R}^{m \times n} $ is easy:
$$ X = \operatorname{Proj}_{\mathcal{T}} \left( Y \right) = \operatorname{tril} \left( Y \right) $$
Namely, zeroing all elements above the main diagonal of $ Y $.
By utilizing the Projected Gradient Descent it is easy to solve this problem:
$$
\begin{align*}
{X}^{k + 1} & = {X}^{k} - \alpha \left( X B {B}^{T} - A {B}^{T} \right) \\
{X}^{k + 2} & = \operatorname{Proj}_{\mathcal{T}} \left( {X}^{k + 1} \right)\\
\end{align*}
$$
The full MATLAB code with CVX validation is available in my StackExchnage Mathematics Q2876283 GitHub Repository.
The solution is very similar to the solution in Q2421545 - Solve Least Squares (Frobenius Norm) Problem with Diagonal Matrix Constraint.
Remark
I think you can also get a closed form solution for each element in $ X $ if you go through deriving the derivative with respect to each element $ X $.
Another approach would be developing the Linear Operator which operates on $ \frac{ \left( n - 1 \right) n }{2} $ elements and creates an $ n \times n $ Triangular Matrix.
My interpretation is that the question is: How can I assure that the linear inequality $q_{ij}-q_{il}-q_{kj}+q_{kl} \geq 0$ holds when $q_{ij} \geq q_{kl}$.
This is an implication betweeen linear inequalities
$$q_{ij} \geq q_{kl} \Rightarrow q_{ij}-q_{il}-q_{kj}+q_{kl} \geq 0 $$
Unfortunately not LP-representable. However, it can be modelled using auxilliary binary variables (thus leading to a mixed-integer quadratic program)
For notational simplicty, consider the generic case for some constraints (not even necessarily linear) $p(x)\geq 0$ and $w(x)\geq 0$ and decision variable $x$
$$p(x) \geq 0 \Rightarrow w(x) \geq 0 $$
Introduce a binary variable $\delta$ and a large number $M$ (more later). The implication is guaranteed to hold if the following two constraints hold
$$p(x) \leq M \delta, w(x)\geq -M(1-\delta)$$
If $p(x)$ is positive, $\delta$ will be forced to be $1$, which forces $w(x)$ to be non-negative.
This is called big-M modelling. When you implement, $M$ should not be made unnecessarily large, but just large enough so that it doesn't affect the feasible space in $x$ when the binary variable is inactive. In your case, assuming the expressions are built from $q$ the difference and sums between your up to 4 $q$ variables can trivially never exceed 4, hence $M=4$ or something like that should be possible. If the variables $a$ and $b$ are involved, you have to know bounds on these so you can bound the expressions accordingly.
If $p(x)$ is a vector of length $n$, and you want the condition to be activated when all elements are non-negative, you would introduce a vector binary variable $\Delta$ in addition to your previous $\delta$, and use $p(x) \leq M\Delta, \delta \geq 1 + \sum\Delta - n$. If all elements are positive, all $\Delta$ binaries will be forced to be $1$, forcing $\delta$ to be $1$.
Best Answer
In case $ B $ is a Positive Definite Matrix then there is $ {C}^{T} C = B $ by the Cholesky Decomposition.
So the problem can be rewritten as:
$$\begin{aligned} \arg \min_{X} \quad & \frac{1}{2} {\left\| A X C \right\|}_{F}^{2} - \operatorname{Tr} \left( D X \right) \\ \text{subject to} \quad & L \leq X \leq U \quad \text{Element wise} \end{aligned}$$
Where $ D = B A $.
Then the gradient of the objective function is easy:
$$ {A}^{T} A X C {C}^{T} - {D}^{T} $$
Now, just use Projected Gradient Descent and you're done.
Remark
In case the matrix is Positive Semi Definite, then you can use the LDL Decomposition and build $ C $ from there in the same manner. If $ B $ is neither PSd nor PD then the problem is not convex. Then you can do the same but only local solution is guaranteed.