\begin{align}
y(t) = \int_{0}^{t} K(t-\tau) u(\tau) d\tau \approx
\int_{0}^{t} \tilde{K}(t-\tau) u(\tau) d\tau =
\begin{cases}
\dot{\mathbf{\xi}}_c(t) = \mathbf{A}_c \mathbf{\xi}_c(t) + \mathbf{B}_c u(t) \\
y(t) = \mathbf{C}_c \mathbf{\xi}_c(t)
\end{cases}
\end{align}

where $K(t)$ is the original impulse or kernel function, $\tilde{K}(t)$ is the approximated impulse function, $\mathbf{\xi}_c$ are the additional states, $\lbrace \mathbf{A}_c,\mathbf{B}_c,\mathbf{C}_c \rbrace$ are the matrices that comprise the state-space that is **exact** with $\tilde{K}(t)$.
The approximate impulse function and the LTI representation are related by the impulse response:
\begin{align}
\tilde{K}(t) = \mathbf{C}_c e^{\mathbf{A}_c t} \mathbf{B}_c
\end{align}

and in the frequency domain:
\begin{align}
\tilde{\mathcal{K}}(s) = \mathbf{C}_c (s\mathbf{I} - \mathbf{A}_c) \mathbf{B}_c = \frac{P_m(s)}{Q_n(s)} = \frac{p_ms^m + p_{m-1} s^{m-1} + \cdots + p_1 s + p_0}{s^n + q_{n-1} s^{n-1} + \cdots + q_1 s + q_0}
\end{align}

where $\tilde{\mathcal{K}}(s)$ is termed the [transfer function](https://en.wikipedia.org/wiki/Transfer_function) (TF) and is a [rational function](https://en.wikipedia.org/wiki/Rational_function). These relationships demonstrate that there is a (nonunique) mapping between a state-space model and transfer function.
---
## Representation
#### Basis Function
Here we will consider the following basis function:
\begin{align}
\phi(t,\mathbf{\theta}) = \theta_1 e^{-\theta_2 t} \cos( \theta_3 t + \theta_4)
\end{align}

which is an exponentially decaying sine wave with phase delay and variable amplitude.
This time-domain function can be represented in the [Laplace](https://en.wikipedia.org/wiki/Laplace_transform) $s$-domain:
\begin{align}
\Phi(s) = \mathcal{L}\{\phi(t)\} =
\theta_1 \frac{\cos(\theta_4) s + \theta_2 \cos(\theta_4) - \theta_3 \sin(\theta_4)}{s^2 + 2\theta_2 s + \theta_2^2 + \theta_3^2} =
\theta_1 \cos(\theta_4) \frac{s+\rho(\mathbf{\theta})}{(s+\zeta(\mathbf{\theta}))(s+\bar{\zeta}(\mathbf{\theta}))}
\end{align}

where the zero and poles of $\Phi(s)$ are:
\begin{align}
-\rho(\mathbf{\theta}) &= -\theta_2 + \theta_3 \tan(\theta_4) \\
-\zeta(\mathbf{\theta}) &= -\theta_2 \pm \theta_3 j
\end{align}

We note that if $\theta_3=0$, then we have a pole/zero cancellation:
\begin{align}
\phi_{\theta_3=0}(t) &= \theta_1 \cos(\theta_4) e^{-\theta_2 t} \\
\Phi_{\theta_3=0}(s) &= \theta_1 \cos(\theta_4) \frac{1}{s+\theta_2}
\end{align}

which is now only an exponential decay with variable amplitude.
#### Function Approximation
The specified impulse response function will be as the sum of $N$ basis functions:
\begin{align}
K(t) \approx \tilde{K}(t) = \sum_{k=1}^N \phi_k(t,\Theta)
\end{align}

where $\Theta$ is the collection of all $\theta$ from each individual basis function. Using the linearity property of the unilateral Laplace transform, the approximation can be written as follows in the $s$-domain:
\begin{align}
\tilde{\mathcal{K}}(s) \approx \tilde{\mathcal{K}}(s) = \sum_{k=1}^N \Phi_k(s,\Theta) = \frac{P_{2N-1}(s,\Theta)}{Q_{2N}(s,\Theta)}
\end{align}

which is a rational function that can be mapped to a state-space model as discussed above.
#### Some Properties of the Chosen Representation
There are some properties of $\tilde{K}(s)$ that should be satisfied to ensure that the approximated system is useful. These properties are now listed and shown to be satisfied by the chosen representation.
- The TF is strictly proper since $m \leq n$. This property is equivalent to:
\begin{align}
\lim_{\omega\to\infty} \tilde{\mathcal{K}}(j\omega) = 0
\end{align}

- The TF's relative degree is 1 since $n - m = 1$.
- The TF is bounded-input, bounded-output (BIBO) stable if all $\theta_{2} > 0$, so we eventually have exponential decay. This property is equivalent to:
\begin{align}
\lim_{t\to \infty} \tilde{K}(t) = 0
\end{align}

- The TF is passive if all the poles have real parts less than zero, i.e., $\theta_{2} > 0$.
---
## Fitting
#### Optimization Problem
Ideally, we want to minimize the square error between the original impulse function and the approximated impulse function over the entire region of interest:
\begin{align}
\min_{\Theta} \ \int_0^{\infty} \left| K(t) - \tilde{K}(t,\Theta) \right|^2 dt = \min_{\Theta} \ \int_0^{\infty} \left| K(t) - \sum_{k=1}^N \phi_k(t,\Theta) \right|^2 dt
\end{align}

However, it is impractical to directly evaluate the optimization problem above so discrete points in time will be evaluated instead to determine the error:
\begin{align}
\min_{\Theta} \ \int_0^{\infty} \left| K(t) - \sum_{k=1}^N \phi_k(t,\Theta) \right|^2 dt
\approx \min_{\Theta} \ \sum_{z=1}^{M} \left| K(t_z) - \sum_{k=1}^N \phi_k(t_z,\Theta) \right|^2
\end{align}

and the selection of the time grid $T$ consisting of $M$ points is important to ensure a suitable final fit.
This fitting problem is the form of the standard nonlinear least-squares problem and MATLAB solvers such as [lsqnonlin](https://www.mathworks.com/help/optim/ug/lsqnonlin.html) can be utilized.
Due to the nonlinear nature of the fitting problem, it can be important to perform a global search procedure to improve the quality of the final solution (such as with [multistart](https://www.mathworks.com/help/gads/multistart.html) or [ga](https://www.mathworks.com/help/gads/ga.html) in MATLAB).
#### Constraints
We also can impose some constraints on the values of $\theta$. Some of the constraints are based on the desired properties discussed above, some are natural constraints, and some are to target the search process to a general region of interest.
\begin{align}
\theta_{1,\min} \leq &\theta_1 \leq \theta_{1,\max} \quad &\text{amplitude} \\
0 < \theta_{2,\min} \leq &\theta_2 \leq \theta_{2,\max} \quad &\text{decay rate} \\
0 \leq &\theta_3 \leq \theta_{3,\max} \quad &\text{frequency} \\
0 \leq &\theta_4 \leq \pi \quad &\text{phase shift} \\
\end{align}

The choice of the bounds is problem specific. The value of $\theta_{2,\min}$ can be selected such that the function decays to (near) zero in a certain amount of time, which is a desired property discussed below. The upper and lower bounds for all the variables are then combined into $L$ and $U$ vectors and are readily included in the nonlinear least-squares solver.
#### Completing the Realization
The partial fraction expansion is combined into the equivalent rational function by incrementally combining each basis function using $P/Q + p/q = (Pq + pQ)/Qq$.
Once the transfer function is determined, we can utilize the [companion form](https://www.mathworks.com/help/control/ref/canon.html#bsq6mpn) directly to create the state-space system.
Alternatively, the [modal form](https://www.mathworks.com/help/control/ref/canon.html#bsq6mph) can be used and is more favorable when performing simulations as it results in smaller state magnitudes.
#### Model Reduction
We can utilize model reduction techniques to potentially reduce the number of states.
This method can be especially beneficial when too large of $N$ is chosen during the fitting process.
States that are removed are based on their Hankel singular values in a balanced realization and have minimal contribution to the value of the output.
The MATLAB functions [modred](https://www.mathworks.com/help/control/ref/modred.html) and [balreal](https://www.mathworks.com/help/control/ref/balreal.html) can be used to perform this task.
Again, we output the state-space model in modal form once the model reduction is completed.
---
## Examples
The sum of three exponentials $0.1 e^{-t/2} + 0.2e^{-2t} + 0.5e^{-3t}$ can be exactly approximated with 3 states:
[![impulse2LTI_ex2](blogs/research/post_2/impulse2LTI_ex2.png)](blogs/research/post_2/impulse2LTI_ex2.png){data-lightbox="blog_imgs"}
The *step-like* impulse function $0.5\tanh(-14(t-0.5))+0.5$ approximated with 6 states:
[![impulse2LTI_ex3](blogs/research/post_2/impulse2LTI_ex3.png)](blogs/research/post_2/impulse2LTI_ex3.png){data-lightbox="blog_imgs"}
A 6 state LTI system representing wave forces:
[![impulse2LTI_ex6](blogs/research/post_2/impulse2LTI_ex6.png)](blogs/research/post_2/impulse2LTI_ex6.png){data-lightbox="blog_imgs"}