perturbed_models.py
# -*- coding: utf-8 -*-
"""
Created on Fri May 24 13:13:53 2013
@author: Simon Streicher
"""
from __future__ import division
import matplotlib.pyplot as plt
import matplotlib as ml
import numpy as np
from Uncertain_MIMO import bound_SISO_wi, weight_calc
from perturbations import possibilities
ml.rcParams['font.family'] = 'serif'
ml.rcParams['font.serif'] = ['Cambria'] + ml.rcParams['font.serif']
"""Calculates pertubated models for specifc elements"""
"""This section is for the g11 element in the process transfer matrix"""
# Define bounds on each uncertain parameter
"""[k, theta, alpha, beta, gamma]"""
umin = [0.8, 13, 0.95, 0.95, 0.95]
umax = [1.2, 17, 1.05, 1.05, 1.05]
#umin = [1, 15, 1, 1, 1]
#umax = [1, 15, 1, 1, 1]
# Define amount of steps to use
steps = 100
# Define start and end of frequency investigated
w_start = -3
w_end = 2
vec = possibilities(umin, umax, 3)
def G(s):
"""Returns value of the nominal plant"""
G = 20.6 * (20 * s + 1) * np.exp(-1 * s * 15) / ((10.0 * s + 1) ** 4)
return G
def Gnp(s):
"""Returns factor part of nominal plant that is not perturbed"""
Gnp = 20.6 / (10.0 * s + 1) ** 2
return Gnp
def Gp(s, vec, i):
"""
Vec contains all perturbations of the varaibles to be used
i is the counter to be used
Gp is the perturbed plant
"""
# Gp = Gnp(s) * ((vec[i, 0] * (vec[i, 4] * 20 * s + 1))
# / ((vec[i, 2] * 10 * s + 1)
# * (vec[i, 3] * 10 * s + 1)))
Gp = Gnp(s) * ((vec[i, 0] * (vec[i, 4] * 20 * s + 1))
/ ((vec[i, 2] * 10 * s + 1)
* (vec[i, 3] * 10 * s + 1))
* np.exp(-1 * s * vec[i, 1]))
return Gp
def weight_i(s, a, b, c):
"""Function to give the weight of the uncertainty"""
w_i = (((s * a + start) / ((a /end) * s + 1))
* ((s ** 2 + b * s + 1) / (s ** 2 + c * s + 1)))
# w_i = ((((s * a + start) / ((a / (end)) * s + 1))
# * ((s ** 2 + b * s + 1) / (s ** 2 + c * s + 1)))
# * ((s ** 2 + d * s + 1) / (s ** 2 + e * s + 1)))
return w_i
# Call bound function
li = bound_SISO_wi(w_start, w_end, vec, G, Gp, steps)
start = np.min(li)
end = np.max(li)
# Call weight calculation function
popt, l_mod, w = weight_calc(w_start, w_end, li, weight_i, steps)
# Modify the weight parameters to get a better fit
popt_mod = popt
popt_mod[0] = popt[0] + 0.5
popt_mod[1] = popt[1]
popt_mod[2] = popt[2]
# Plot the weight and the li error it must be above
# The red line must be above the blue line at all time
plt.loglog(w, l_mod, 'b')
plt.loglog(w, weight_i(w, popt_mod[0], popt_mod[1], popt_mod[2]), 'r')
plt.show()