#!/usr/bin/env python
# coding: utf-8

# In[1]:


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


# In[2]:


get_ipython().run_line_magic('matplotlib', 'nbagg')
# only needed when developing python-control
get_ipython().run_line_magic('load_ext', 'autoreload')
get_ipython().run_line_magic('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)


# ## 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.))


# 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


# ## 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]