The first reference:
 System Identification — Theory For the User, Lennart Ljung, Section 13.4-13.5, 2nd ed, PTR Prentice Hall, Upper Saddle River, N.J., 1999.
explains the identification in closed loop quite well. That is the basic perspective taken in System Identification Toolbox. To answer the questions raised by OP:
Should I use a MISO or MIMO approach to system identification?
It does not matter. As in the case of open-loop identification, it is sometimes advisable to split a multi-output estimation into several single output-estimations. This reduces the number of parameters seen by the nonlinear least squares optimizer and therefore reduces the chances of getting stuck in a local minima. The main thing with closed-loop id is to use a structure that has a rich noise structure; ssest , bj, armax commands are some options.
How can I best identify the unstable plant and verify the model using the system identification toolbox?
Unstable id is fundamentally hard since any small errors owing to noise on wrong structure are amplified. Direct estimation of unstable plant using time-domain data rarely works well. Some tips:
- If you must use time-domain data, use focus = 'prediction' option (which is the default anyway) on a model structure with non-trivial noise components (e.g., SS, BJ, ARMAX). It can happen sometimes that the estimation is able to find a stable predictor for an unstable model (abs(eig(A))>1 but abs(eig(A-K*C))<1 in a state-space realization). If this happens you have found your unstable model by time-domain identification.
- Another approach is to perform sine-wave experiments in closed-loop (measuring u and y as before) and use it to build a frequency response of the model. If the original model is in Simulink, you can use the frestimate command of Simulink Control Design Toolbox (with sine-stream inputs). Then use the estimated frequency response (packaged as FRD object) to perform a frequency-domain identification. In frequency-domain, identification of unstable models works out significantly more easily. Somehting like:
g = frestimate(model, inputs,...
m = ssest(g, nx);
When the plant is unstable, how do I simulate the model and compare to verification data sets?
Unstable simulations is inadvisable. Your best bet is to compare responses in frequency-domain. Either use the FRD object or use FFT(time_domain_data) as reference data in COMPARE command.