# Numerical computations for transfer functions / state space models?

Let's say that I have the transfer function:

$$G(s) = \frac{1}{s^2 + s + 1}$$

And I want to convert that to a state space model, view the frequency response in a bode diagram, plot the transfer function in a nyquist diagram and check the stability etc...

Or If I want to transfer a MIMO state space model to a transfer function matrix?

I know how to do that by hand. But is there a way to do that numerically by using pure clean MATLAB / Octave code? No Toolbox allowed. Only built in functions.

Let's say I going to compute the state transition matrix :

$$\Phi(t) = L^{-1}(sI-A)^{-1}$$

I know there is a numerical way:

$$\Phi = I + At + \frac{A^2t^2}{2!}+ \frac{A^3t^3}{3!}+ \frac{A^4t^4}{4!} \dots$$

Or compute the transfer function from a state space model:

$$G(s) = C(sI-A)^{-1}B + D$$

Or TF -> SS...I don't remember any numerical method for that. I know for sure that everything I do will be a state space model. If I going to merge two transfer functions, I need to transform both of them into SS and then merge them together, then SS -> TF.

If I want to check the step answer for a transfer function, I need to convert that to TF -> SS and then simulate. But how?

Edit: Should I change language? Julia? Fortran 77?

## 1 Answer

Writing a SISO transfer function as a state space model can be done as follows. Given the following general SISO transfer function

$$G(s) = \frac{a_0 + a_1\,s + \cdots + a_{n-1}\,s^{n-1} + a_n\,s^n}{b_0 + b_1\,s + \cdots + b_{n-1}\,s^{n-1} + s^n}$$

an equivalent state space model in the controllable canonical form can be found with

$$\left[\begin{array}{c|c} A & B \\ \hline C & D \end{array}\right] = \left[\begin{array}{ccccc|c} 0 & 1 & 0 & \cdots & 0 & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 & \cdots & 0 & 1 & 0 & \vdots \\ 0 & \cdots & \cdots & 0 & 1 & 0 \\ -b_0 & -b_1 & \cdots & -b_{n-2} & -b_{n-1} & 1 \\ \hline a_0-a_n\,b_0 & a_1-a_n\,b_1 & \cdots & a_{n-2}-a_n\,b_{n-2} & a_{n-1}-a_n\,b_{n-1} & a_n \end{array}\right]$$

In order to obtain a Bode or Nyquist diagram you just have to evaluate $G(j\,\omega)$ for a range of value for (angular) frequencies $\omega$. For the Bode diagram you have to calculate the dB of the magnitude and the argument of this result and plot these against $\omega$ on a logarithmic scale. For the Nyquist diagram you can directly plot the real and imaginary part of that result. For stability you can just look at the roots of the characteristic polynomial or the eigenvalues of the system matrix.

For example the following code could be used:

N = 5;                      % System order
a = rand(1, N);             % Numerator coefficients
b = rand(1, N); b(N) = 1;   % Denominator coeffcients

L = 1e3;                    % Number of frequency elements
w = logspace(-3, 3, L);     % Angular frequencies

% Evaluate transfer function
G = zeros(1, L);
for k = 1 : L
G(k) = (a * (1i * w(k)).^(0 : N-1).') / (b * (1i * w(k)).^(0 : N-1).');
end

% Plot Bode diagram
figure
subplot(2,1,1)
semilogx(w, 20*log10(abs(G)))
subplot(2,1,2)
semilogx(w, angle(G) * 180/pi)

% Plot Nyquist diagram
figure
plot([real(G) nan real(G)], [imag(G) nan -imag(G)], -1, 0, '+')

• But what if I don't know the state space form and the transfer funtion look like a mess? – Daniel Mårtensson Sep 5 '17 at 16:48
• What happen if the the TF model look like this: $$G(s) = \frac{a_0 + a_1\,s + \cdots + a_{n-1}\,s^{n-1} + a_n\,s^n}{b_0 + b_1\,s + \cdots + b_{n-1}\,s^{n-1} + s^n}e^{-sT}$$ – Daniel Mårtensson Sep 5 '17 at 20:40
• Yes I know! It will be a discrete time state space model. But the issue is...how can I make so MATLAB/Octave can count the numerators and denumerators if I just type in a random linear TF? I don't asking you to solve the problem for me. I just asking for guidance. – Daniel Mårtensson Sep 5 '17 at 20:54
• @DanielMårtensson If your transfer function is in the s domain but also contains delays (or other not mentioned functions) you can still evacuate it numerically at various frequencies in order to plot the Bode and Nyquist diagram. For closed loop stability analysis you can still use the Nyquist diagram, but there will be infinitely many poles. When you want to convert it into a state space model you either would need to discretize (I believe the discretization does not need to have the same sampling time as the delay). Or you can use a Padé approximation. – Kwin van der Veen Sep 5 '17 at 22:08
• Thanks. If I use the command tf and set s = tf ('s') . The problem for men is that I need to create a symbolic variable s in Octave/MATLAB if I don't have any control package. How can I tell octave that s is a laplace s and if I with this "G = 3/(s^2 + s)". Octave are going to recognize and count how many poles and zeros they are. – Daniel Mårtensson Sep 6 '17 at 5:39