Welcome to Robotics, PaoloH! This is a fantastic question for Robotics - It has some Matlab/Simulink, some control theory, some spatial (quaternion) representations, etc. Robotics is the place to come when your question spans multiple fields!
In looking at your question, the thing that I noticed is that your reference quaternion is $[0; 1; 0; 1]$. It is not a unit quaternion, and I believe this may be your issue.
I looked over your block diagram, and I didn't see anything glaringly wrong there. As SteveO mentioned, the way you're treating the references is a little unusual, but the math all works out. I can't see what you're doing behind the "reference quaternion" or "error quaternion" blocks, but let's take a look at that unit quaternion.
Right now, the magnitude of your reference quaternion is
$$
\sqrt{0^2 + 1^2 + 0^2 + 1^2} = \sqrt(2) \approx 1.414 \\
$$
If you want to convert your reference quaternion to a unit quaternion, then you divide each term of the quaternion by the magnitude of the quaternion, and you wind up with a reference unit quaternion of:
$$
q_{\mbox{ref}} = \left[\begin{array}{ccc}
0 \\
\frac{1}{\sqrt{2}} \\
0 \\
\frac{1}{\sqrt{2}} \\
\end{array}\right] \\
q_{\mbox{ref}} \approx \left[\begin{array}{ccc}
0 \\
0.707 \\
0 \\
0.707 \\
\end{array}\right];
$$
You can review your quaternion output plot and see that q(2) and q(4) are both moving toward a numeric value of ~0.7. The only real problem seems to be that the polarity on q(2) is wrong.
I would guess your "reference quaternion" block is to make the skew symmetric matrix for the "error quaternion" block? The sign problem on your quaternion output could be hiding anywhere, but I'd double check the skew symmetric matrix (that it is actually skew symmetric; $\Omega^T == -\Omega$), and then I'd check the gains. When I said the math all works out on the reference handling, I double checked that just now and it works out for the $\omega$ speed handling. I can't actually tell for the quaternion handling.
Typically, error is (reference - feedback), then you apply gains and sum the scaled error terms.
You have your error as (feedback - reference), and then you apply gains and negate the scaled error terms. BUT, it looks like, for quaternion error, you actually ARE taking (reference - feedback) but then you're still inverting it anyways.
If this isn't enough to get your question resolved, please edit your question to show what's going on under the "Reference quaternion" and "error quaternion" blocks.