R isn’t really the first programming language engineers learn or choose to use. Matlab is the de facto standard for everything involving controller design. But not all control theoretical tools are solely used by control engineers. Frequency-domain response analysis (FdRA) for example is used by systems biologists to either identify the structure of the model from inputs and outputs alone or to tune output behavior by changing the input. Thus, I set to implement basic FdRA functions in R and recently published as freqr.

Here I shortly wanted to write some background information on FdRA to go with the R package.

What is FdRA?

Any linear time-invariant system with multiple inputs and multiple outputs (MIMO) can be given in state-space representation such that

with \(\boldsymbol{A}\in\mathbb{R}^{n\times n}\), \(\boldsymbol{B}\in\mathbb{R}^{n\times p}\), \(\boldsymbol{C}\in\mathbb{R}^{q\times n}\), and \(\boldsymbol{D}\in\mathbb{R}^{q\times p}\). For a system with a single input (\(p=1\)) and a single output \(q=1\) (SISO) this reduces to

The transfer function now follows with \(\boldsymbol{x}(0)=\boldsymbol{x}_0\) to

with \(\boldsymbol{G}\in\mathbb{R}^{q\times p}\). In the SISO case the transfer function reduces to \(G\in\mathbb{R}\). Let’s first just consider the SISO case.

SISO systems

If a LTI system with a transfer function \(G(s)\) is excited by the harmonic input

then following attenuation of all transient processes the output \(y(t)\) of the system will again be a harmonic function with the identical frequency \(\omega\) but a different phase and amplitude

Bode plot

The input/output behavior of a system represented by its transfer function can now be visualized with the help of a Bode plot. It is composed of

1. the magnitude in form of the graph \((\log(\omega), \lvert G(i\omega)\rvert_\text{dB})\) and
2. the phase in form of the graph \((\log(\omega), \arg(G(i\omega)))\) where the phase is given in degree.

Note that the magnitude in decibels can be calculated from

The magnitude plot now informs on whether the input amplitude is amplified or attenuated. The phase plot on the other can provides information on the time scale the system operates on, i.e. how fast it reacts to changes in the input.

MIMO systems

For a system with multiple inputs and multiple outputs obtaining the frequency response is not so easy. Because the inputs and outputs are now vectors instead of scalars, their direction matters as well. The singular value decomposition (SVD) provides a useful way of quantifying this multivariable directionality [¹]. Note, that in the following \(u\) will not denote the inputs but the left-singular vectors of the singular value decomposition.

So, the transfer matrix \(\boldsymbol{G}(s)\) may now be decomposed into its singular values as

where the orthonormal column vectors of \(\boldsymbol{U}\), denoted as \(\boldsymbol{u}_i\), represent the output directions of the system while the orthonormal column vectors of \(\boldsymbol{V}\), denoted as \(\boldsymbol{v}_i\), represent the input directions.

Since \(\boldsymbol{V}\) is unitary it follows that for column \(i\) we can write

which means that for an input in direction \(\boldsymbol{v}_i\) the output is in direction \(\boldsymbol{u}_i\).

The maximum gain for any input direction is equal to the maximum singular value

while the smallest gain for any input direction is equal to the minimum singular value

with \(k=\min(p,q)\).

If we now define \(\boldsymbol{u}_1=\overline{\boldsymbol{u}}\), \(\boldsymbol{v}_1=\overline{\boldsymbol{v}}\), \(\boldsymbol{u}_k=\underline{\boldsymbol{u}}\), and \(\boldsymbol{v}_k=\underline{\boldsymbol{v}}\), it follows that \(\boldsymbol{G}\overline{\boldsymbol{v}}=\overline{\sigma}\,\overline{\boldsymbol{u}}\) and \(\boldsymbol{G}\underline{\boldsymbol{v}}=\underline{\sigma}\,\underline{\boldsymbol{u}}\). The vector \(\overline{\boldsymbol{v}}\) corresponds to the input direction with the largest amplification, and \(\overline{\boldsymbol{u}}\) is the corresponding output direction in which the inputs are most effective.

Sigma plot

The singular values and with them the maximum and minimum gain can now be visualized in a singular value plot.

freqr

freqr now allow you to generate exactly those two important visualization of the frequency response by calling either bodeplot or sigmaplot.

To install freqr use

devtools::install_github("pascalschulthess/freqr")

For more details on the package and on usage and examples please visit the corresponding GitHub repository.

References