Numerical computations for transfer functions / state space models?

Question

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?

Answer 1

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, '+')