I was given a state-space model of a system and a list of outputs for t=0 to t=5, sampled every 0.1 seconds and asked to approximate the initial condition. Is there a function in MATLAB that will take a long list of data and approximate the initial condition of the system?
The system $$\dot x=\begin{bmatrix}-1 & 1 & 0\\0 & -1 & -1\\1 & 0 & -2\end{bmatrix}x + \begin{bmatrix} 0 & 1\\1 & 0\\0 & 0\end{bmatrix}u$$
$$y = \begin{bmatrix}1 & 0 & 1\end{bmatrix}x$$
was given the input $u(t) = \begin{bmatrix}0\\0\end{bmatrix}$ for $t = 0$ to $5$ seconds and the following data was obtained by sampling every 0.1 seconds. Approximately what was the initial condition $x(0)$?
$$\begin{array}{c|l} t & \text{y(t)}\\ \hline 0.0 & 0.0000\\ 0.1 & 0.556011\\ 0.2 & 1.02816\\ 0.3 & 1.42273\\ 0.4 & 1.74606\\ 0.5 & 2.00443\\ 0.6 & 2.20402\\ 0.7 & 2.35082\\ 0.8 & 2.45059\\ 0.9 & 2.50882\\ 1.0 & 2.5307\\ 1.1 & 2.52112\\ 1.2 & 2.48461\\ 1.3 & 2.4254\\ 1.4 & 2.34736\\ 1.5 & 2.25403\\ 1.6 & 2.14863\\ 1.7 & 2.03404\\ 1.8 & 1.91286\\ 1.9 & 1.78739\\ 2.0 & 1.65964\\ 2.1 & 1.53138\\ 2.2 & 1.40411\\ 2.3 & 1.27914\\ 2.4 & 1.15753\\ 2.5 & 1.04019\\ 2.6 & 0.927817\\ 2.7 & 0.82098\\ 2.8 & 0.720091\\ 2.9 & 0.625442\\ 3.0 & 0.53721\\ 3.1 & 0.455476\\ 3.2 & 0.380238\\ 3.3 & 0.311418\\ 3.4 & 0.248879\\ 3.5 & 0.192432\\ 3.6 & 0.141847\\ 3.7 & 0.0968592\\ 3.8 & 0.0571798\\ 3.9 & 0.0225006\\ 4.0 & -0.00749888\\ 4.1 & -0.0331464\\ 4.2 & -0.0547719\\ 4.3 & -0.0727038\\ 4.4 & -0.0872657\\ 4.5 & -0.0987728\\ 4.6 & -0.10753\\ 4.7 & -0.11383\\ 4.8 & -0.11795\\ 4.9 & -0.120156\\ 5.0 & -0.120694\\ \end{array} $$
Let's make this generic: consider the linear system $\dot{x}=Ax+Bu$, $y=Cx+Du$. If $u\equiv 0$ then the specific values of $B,D$ are irrelevant and $$x(t)=e^{At}x_0, \quad y(t) = Cx(t) =Ce^{At}x_0$$ where $x_0$ is the initial state. In your case, you have data for $$t\in\{0,T,2T,\dots,NT\} \qquad T=0.1, \quad N=50$$ If we define $\bar{A}\triangleq e^{AT}$ and $y_k\triangleq y(kT)$, then $y_k = CA^kx_0$. Therefore $$\begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_N \end{bmatrix} = \begin{bmatrix} C \\ C\bar{A} \\ C\bar{A}^2 \\ \vdots \\ C\bar{A}^N \end{bmatrix} x_0.$$ Construct this overdetermined $(N+1)\times n$ linear system and solve with a simple backslash operation to find $x_0$. If the measurements are exact and there is no noise in the system, the value of $x_0$ obtained will be exact as well.
In practice, neither the assumption of exact measurements or noise-free operation is valid. Therefore, I'd weight the earlier samples more than the later ones. Assuming no system noise and exact measurements, this would give the same result for any $w_0,w_1,\dots,w_N>0$: $$\begin{bmatrix} w_0 y_0 \\ w_1 y_1 \\ w_2 y_2 \\ \vdots \\ w_N y_N \end{bmatrix} = \begin{bmatrix} w_0 C \\ w_1 C\bar{A} \\ w_2 C\bar{A}^2 \\ \vdots \\ w_N C\bar{A}^N \end{bmatrix} x_0.$$ A common choice is $w_k=\alpha^k$ for some forgetting factor $0<\alpha<1$. A more rigorous development would use statistics about the noise and/or measurement errors to come up with proper weighting. Indeed, that's the very kind of rigor you find in Kalman filter theory.