In [1]:

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import control as ct

In [2]:

%matplotlib nbagg
# only needed when developing python-control
%load_ext autoreload
%autoreload 2

Define continuous system¶

In [3]:

# Distillation column model as in Equation (3.81) of Multivariable Feedback Control, Skogestad and Postlethwaite, 2st Edition.

den = [75, 1]
G = ct.tf([[[87.8], [-86.4]],
          [[108.2], [-109.6]]],
         [[den, den],
          [den, den]])
display(G)

$$\begin{bmatrix}\frac{87.8}{75 s + 1}&\frac{-86.4}{75 s + 1}\\\frac{108.2}{75 s + 1}&\frac{-109.6}{75 s + 1}\\ \end{bmatrix}$$

Define sampled system¶

In [4]:

sampleTime = 10
display('Nyquist frequency: {:.4f} Hz, {:.4f} rad/sec'.format(1./sampleTime /2., 2*np.pi*1./sampleTime /2.))

'Nyquist frequency: 0.0500 Hz, 0.3142 rad/sec'

In [5]:

# MIMO discretization not implemented yet...

Gd11 = ct.sample_system(G[0, 0], sampleTime, 'tustin')
Gd12 = ct.sample_system(G[0, 1], sampleTime, 'tustin')
Gd21 = ct.sample_system(G[1, 0], sampleTime, 'tustin')
Gd22 = ct.sample_system(G[1, 1], sampleTime, 'tustin')

Gd = ct.tf([[Gd11.num[0][0], Gd12.num[0][0]],
            [Gd21.num[0][0], Gd22.num[0][0]]],
           [[Gd11.den[0][0], Gd12.den[0][0]],
            [Gd21.den[0][0], Gd22.den[0][0]]], dt=Gd11.dt)
Gd

Out[5]:

$$\begin{bmatrix}\frac{5.487 z + 5.488}{z - 0.875}&\frac{-5.4 z - 5.4}{z - 0.875}\\\frac{6.763 z + 6.763}{z - 0.875}&\frac{-6.85 z - 6.85}{z - 0.875}\\ \end{bmatrix}\quad dt = 10$$

Draw Singular values plots¶

Continuous-time system¶

In [6]:

omega = np.logspace(-4, 1, 1000)
plt.figure()
response = ct.freqplot.singular_values_response(G, omega)
sigma_ct, omega_ct = response
response.plot();

Discrete-time system¶

In [7]:

plt.figure()
response = ct.freqplot.singular_values_response(Gd, omega)
sigma_dt, omega_dt = response
response.plot();

Continuous-time and discrete-time systems altogether¶

In [8]:

plt.figure()
ct.freqplot.singular_values_plot([G, Gd], omega);

Superposition on the same singular values plot¶

In [9]:

plt.figure()
ct.freqplot.singular_values_plot(G, omega);

In [10]:

ct.freqplot.singular_values_plot(Gd, omega);

Analysis in DC¶

In [11]:

G_dc = np.array([[87.8, -86.4],
                 [108.2, -109.6]])

In [12]:

U, S, V = np.linalg.svd(G_dc)

In [13]:

S, sigma_ct[:, 0]

Out[13]:

(array([197.20868123,   1.39141948]), array([197.20313497,   1.39138034]))