Obviously, we had better make \(X\) non-decreasing. In this editorial, we implicitly assume that \(X\) is non-decreasing.
First, we may assume that we immediately discard the integers that does not end up in \(X\). Under this assumption, \(S\) should end up in being empty. Moreover, any integer pushed to \(S\) should be not less than any other elements in \(S\); especially, when the maximum integer \(M\) in the resulting \(X\) is moved to \(S\) for the first time, \(S\) should be empty.
Based on this idea, consider the following \(\mathrm {dp}\). (DP stands for Dynamic Programming.)

• \(\mathrm {dp}_{i,j,k,l}\) : the maximum score when obtaining \(X\) using \(a_i,a_{i+1},\ldots,a_j\) that are between \(k\) and \(l\) (inclusive)

The answer for this problem is \(\mathrm{dp}_{1,N,1,50}\).
Now we consider the transitions.

Transition when \(l\) is not contained in resulting \(X\)

It corresponds to \(\mathrm {dp}_{i,j,k,l-1}\).

Transition when \(l\) is contained in resulting \(X\)

For each term such that \(a_m=l\), simulate the transition assuming that \(l\) is moved for the first time when moving that term. As mentioned, \(S\) must be empty at this time. Also, if \(n\) is the maximum term that occurs before \(a_m\) in the resulting \(X\), then we cannot include an integer less than \(n\) after \(a_m\) in \(X\). (Note that \(n\) is already the last term of \(X\) because \(S\) is empty.)
Conversely, let \(X_1\) be the part of \(X\) that comes before \(a_m\), and \(X_2\) be the part that comes after \(a_m\) (without using integers less than \(n\)); then we can obtain \(X\) of \((1+|X_1|+|X_2|)\) terms by concatenating \(a_m,X_1,X_2\). Therefore, the transition \(\mathrm{dp}_{i,j,k,l} = \underset{k \leq n \lt l}{\max}(1 + \mathrm{dp}_{i,m-1,k,n} + \mathrm{dp}_{m+1,j,n,l})\) is valid. (Due to the discussion so far in which we assumed that \(m\) is the first position such that \(a_m=l\), we don’t include \(l\) in the range of \(n\), but including it is OK)
We analyze the time complexity. We assume that the maximum value of \(a_i\) equals \(N\). (Even if not, we can do coordinate compression to make it \(N\).) There are \(\mathrm{O}(N^4)\) states of \(\mathrm{dp}\), and in a naive implementation, we need to find the positions such that \(a_m=l\) for each DP element, so this part costs \(\mathrm{O}(N^5)\) time. Moreover, inspecting \(n\) for each \(a_m=l\) imposes an \(\mathrm{O}(N)\) transition, which occurs when \(l=a_i,l=a_{i+1},\ldots,l=a_{j}\) for each \((i,j,k)\) so it costs a total of \(\mathrm{O}(N^5)\) time too. Directly assigning the constraints into this polynomial yields \(50^5=312500000\), which looks too large for a time limit of two seconds, but the constant factor is small enough to meet the time limit. Some implementation may require \(\mathrm{O}(N^6)\), which may still fit in the time limit if you are careful enough of the constant factors.