Some points:
Don't do "algebra" with tf variables, use the functions series, parallel and feedback.
sys_mimo1 = ss(sys_mimo) is going to convert the tf model to a ss model, I don't think there is a need for the extra step with ssdata.
This one here
sys_mimo_tf(1,2) = sys_mimo_tf(1,2)*feedback(sys_mimo_tf(1,2),C);
is probably not doing what you want it to, it will give you the transfer function
sys_mimo_tf(1,2)^2 / 1 + sys_mimo_tf(1,2)*C
so that's probably where many extra states come from. You should expect 2 extra states per added PID.
Lastly, look at the Hankel Singular Values or a similar measure to figure out whether there might be a lot of almost unobservable/almost uncontrollable states in your resulting model. That would indicate numerical issues.
Best Answer