How do I interpret states in the MATLAB System Identification Toolbox?

Due to the fact that you are mentioning the $A,B,C,D$ matrices I assume that your model is a linear one of the form:

$$ dot{x} = Ax + Bu$$ $$ y = Cx + Du $$

where $x$ is the state vector and $y$ the output of the system. One thing you can do, at least I would go like this, is transform your system to the $text{s-domain}$, that is obtain the input to output transfer function:

$$ T(s) = frac{U(s)}{Y(s)} $$ Since you can measure both the input and the output of the system, collect some data by stimulating the system with some input signal and measure the output. Store both input and output data to vectors. The input signal though should not be some random arbitrary signal. It should be a collection of step inputs, a pulse signal, a summation of sinusoids or even better a random sequence such as (pseudorandom binary sequence or gaussian white noise). The sum of sinusoids can excite a finite number of frequencies, particularly the signal:

$$ u = sum_{i=1}^m A_isin(ω_it) $$

is called sufficiently rich of order $n$ ($n$ is the number of the system's unknown parameters) if and only if $m geq frac{n}{2}$. So, in order to excite all the necessary frequencies of your system you need to know how many unknown parameters you need to estimate. On the other hand, random sequences as the ones I referred to previously, can indeed excite all the frequencies and could be a really good option in order to obtain input-output data which contain all the information needed to perform a really good system identification. I really encourage you to take a look at the Robust Adaptive Control book, which is available online. It really deals with those issues in a great in-depth way.

Now, you have the equation of the transfer function (it should be a third order) and a collection of input and output data. You got everything you need in order to proceed and perform the identification of the unknown parameters of the system.

I find it easier to use the System Identification Toolbox in order to estimate transfer function models. When estimating transfer function models you only need to specifically set the input and output of the system. Of course do not forget the very important step of processing your data. You may have to attenuate any possible noise, choose which part of the data you would like to use (because you have a linear system it would be best to work with output data that present a linear behavior with respect to the input data), normalize them or a whole bunch of them. Of course the toolbox can do these stuff for you by selecting the appropriate choices.

To conclude, you have now obtained a transfer function model of your system. You can easily convert it into a state-space representation and apply any state feedback control law you want.

When identifying a system, you effectively compute a system that happens to have the same response to your inputs as the physical system. However, this does not mean that the states do physically represent the same. You find a transformed system that can be represented as such: $$Tdot{x} = (TAT^{-1})Tx + TBu$$ $$y = (CT^{-1})Tx + Du$$ where $T$ is a transformation matrix. The eigenvalues of $TAT^{-1}$ are equal to the eigenvalues of $A$. And if im correctly, the corresponding transfer functions are equal as well. Retrieving this transformation matrix $T$ is impossible without knowing the actual system. Luckily, you dont have to! As long as you use the same inputs and outputs when identifying the system and when you control it, the controller will work. If you have more inputs or outputs, you must make sure the vector does not get scrambled (so $u_1$ stays $u_1$). If you want to compute a state feedback gain $K$ using these transformed model matrices, you will get $KT^{-1}$ instead. This is required because the state observer will not observe $x$ but $Tx$.

EDIT:

Due to the size of my answer, I'd like to answer your questions here. First of all, I strongly advice against estimating states by either manually differentiating or integrating measurable states. If you differentiate your measured signal, you will also differentiate sensor noise. Due to the unpredictable (yet small) changes in noise, the differentiated signal will amplify the sensor noise. Integrating the sensor bias will ensure the result will drift away from the actual value. To obtain an estimate of states you cannot measure directly, you should use a State observer. Presumably the most popular observer is the Luenberger observer and it works as follows: suppose you create a state estimate called $hat{x}$, which has the following state update: $$dot{hat{x}} = Ahat{x} + Bu + L(y - Chat{x})$$ Where $A, B, C$ are the estimated model matrices, $u$ is the input you put in the system and $y$ is the measured output. The estimation error can be represented as: $$e = hat{x} - x, quad dot{e} = dot{hat{x}}- dot{x}$$ By substituting the equations for $dot{hat{x}}$ and $dot{x}$, the estimation error can be rewritten to: $$dot{e} = (A-LC)e$$ To ensure the error will decrease to 0, one has to create $L$ such that the eigenvalues of $A-LC$ are less than 0. In practical sense, $L$ is chosen such that $A-LC$ converges much faster than $A-BK$. So how does one implement this? Usually one likes to avoid differentiating manually inside code so you would like to convert this system to a discrete system. Matlab commands like c2d.m can do that in a very lazy manner. Next, you implement the observer into arduino using a dummy state $hat{x}$ with 0 as initial value. For each iteration, you compute the update of $hat{x}_{k+1}$ and use the current value of $hat{x}_k$ to compute the next input ($u= -Khat{x}_k$). Feed this input to the system, measure its response and repeat.