How to minimize a weighted sum of RMSE-like terms?

I am trying to solve the following problem: \begin{align} \min&\quad f(x) = \sum_{i=1}^{n}{a_ix_i} + \sum_{i=1}^{n}{b_i\sqrt{\sum_{j=1}^{m}{\left(y_{i,j}-x_i\right)^2}}}\\\text{s.t.}&\quad x_{i+1}\geq x_i \quad\forall i \in \{1,\dots,n-1\} \\&\quad y_{i,j} \geq x_i \quad\forall (i,j) \in \{1,\dots,n\} \times \{1,\dots,m\} \\&\quad y_{i,j} \geq c_j \quad\forall (i,j) \in \{1,\dots,n\} \times \{1,\dots,m\} \\&\quad x_i \geq 0 \quad\forall i \in \{1,\dots,n\} \end{align}

The parameters $$a_i$$, $$b_i$$ and $$c_j$$ are all positive numbers.

The problem is not quadratic because of the $$\sqrt{\cdot}$$ functions in the objective. Therefore I believe mathematical programming solvers will be unable to solve this problem. Are local search methods the way to go here or is there a way to reformulate my problem that makes it easier to solve? If it's the former, are there off-the-shelf tools that I can use ?

The square root terms are two-norms and the rest of the problem is linear. Therefore, this can be solved by a solver capable of handling Second Order Cone Problems (SOCP). Among such solvers are CPLEX, Gurobi, Mosek, SeDuMi, SDPT3, and KNITRO, among others.

This problem is convex, and can be solved with such solvers to global optimality.

This problem can be expressed as a SOCP by using an epigraph formulation of norm. That is, replace the norm terms by $$t_i$$, and introduce constraints, $$\text{norm}(...) \le t_i$$.

If a convex optimization modeling tool such as CVX, YALMIP, CVXPY, or CVXR is used, this reformulation can be handled automatically (under the hood) by the tool.

In CVX, the formulation would look like

cvx_begin
variables x(n) y(n,m)
minimize(a'*x + b'*norms(y - repmat(x,1,m)),2,2)
<incorporation of linear constraints left as straightforward exercise>
cvx_end

making use of norms, which is a CVX function which computes a vector of norms at once from a matrix and produces a vector output. But this can be done using loops using plain old norm instead.

help norms

norms Computation of multiple vector norms. norms( X ) provides a means to compute the norms of multiple vectors packed into a matrix or N-D array. This is useful for performing max-of-norms or sum-of-norms calculations.

All of the vector norms, including the false "-inf" norm, supported
by NORM() have been implemented in the norms() command.
norms(X,P)           = sum(abs(X).^P).^(1/P)
norms(X)             = norms(X,2).
norms(X,inf)         = max(abs(X)).
norms(X,-inf)        = min(abs(X)).
If X is a vector, these computations are completely identical to
their NORM equivalents. If X is a matrix, a row vector is returned
of the norms of each column of X. If X is an N-D matrix, the norms
are computed along the first non-singleton dimension.

norms( X, [], DIM ) or norms( X, 2, DIM ) computes Euclidean norms
along the dimension DIM. norms( X, P, DIM ) computes its norms
along the dimension DIM.

Disciplined convex programming information:
norms is convex, except when P<1, so an error will result if these
non-convex "norms" are used within CVX expressions. norms is
nonmonotonic, so its input must be affine.