In [ ]:
import cvxopt as cvx
import control as con

import numpy as np

In [ ]:

Parameters of some Examples from different sources


In [ ]:
##############################
# Boris Paper UBoot          #
##############################
A = np.matrix([[0., 1., 0.],
              [0., 0., 1.],
              [0., 0., -0.005]])

############### TODO ###################
############### CHECK ##################
#last line of Regelungsnormalform (characteristical polynom coefficients) in a colunm vector
a = -np.matrix(A[np.array(A).shape[1]-1])[np.newaxis].T 

b = np.matrix([[0], [0], [1.]])

d = 0
c = np.matrix([[1], [0], [0]]).T

print np.matrix(con.matlab.pole(con.matlab.ss(A, b, c, d))).T

u_max = 2.5e-5

X0 = [np.matrix([-10, -0.05, -0.0046]).T,
      np.matrix([-10, -0.05,  0.0046]).T,
      np.matrix([-10,  0.05, -0.0046]).T,
      np.matrix([-10,  0.05,  0.0046]).T,
      np.matrix([ 10, -0.05, -0.0046]).T,
      np.matrix([ 10, -0.05,  0.0046]).T,
      np.matrix([ 10,  0.05, -0.0046]).T,
      np.matrix([ 10,  0.05,  0.0046]).T]

x_max = [10, 0.05, 0.0046]

#simulation time
T = 1.0/200.0

#Introduced for fun
u_max_sys = u_max

a_hat      = np.matrix([[4.4469e-8], [2.3073e-5], [4.9148e-3]])
a_hat_star = np.matrix([[1.073e-7], [4.919e-5], [10.4078e-3]])

R1 = np.matrix([[1.6021e-5, 3.26098e-3, 0.4031],
               [3.2698e-3, 1.5666,     163.46],
               [0.4031,    163.46,     40.713]])

##### Entwurf parameter #####
beta = 2 # beta >=1 !

In [ ]:
##############################
# Boris Diss Fusionsreactor  #
##############################
A = np.matrix([[1, 0],
              [0, -0.5]])
b = np.matrix([[1],[-0.5]])
c = np.matrix([1, 0])
u_max = 1
n = b.size

print "Poles:\n", np.matrix(con.matlab.pole(con.matlab.ss(A, b, c, 0))).T

X0 = [np.matrix([-0.9, -2.8]).T,
      np.matrix([-0.9, 2.8]).T,
      np.matrix([0.9, -2.8]).T,
      np.matrix([0.9, 2.8]).T]

x_max = [0.9, 2.8]

#print "A:\n", A
#print "a:\n", a
#print "b:\n", b
#print "c:\n", c

##### Entwurf parameter #####
beta = 2 # beta >=1 !

In [ ]:
###################################################
# Boris Diss Hydraulisches Positioniersystem S.61 #
###################################################
A = np.matrix([[0,      1,        0],
               [0, -0.126,     83.5],
               [0,  -6.89, -0.00707]])
 
b = np.matrix([[0],[0],[-55.51]])
c = np.matrix([1, 0, 0])

print "Poles:\n", np.matrix(con.matlab.pole(con.matlab.ss(A, b, c, 0))).T

u_max = 1
n = len(b)

X0 = [np.matrix([-13.38, -10.7, -2.58]).T,
      np.matrix([-13.38, -10.7,  2.58]).T,
      np.matrix([-13.38,  10.7, -2.58]).T,
      np.matrix([-13.38,  10.7,  2.58]).T,
      np.matrix([ 13.38, -10.7, -2.58]).T,
      np.matrix([ 13.38, -10.7,  2.58]).T,
      np.matrix([ 13.38,  10.7, -2.58]).T,
      np.matrix([ 13.38,  10.7,  2.58]).T]

x_max = [13.38, 10.7, 2.58]

#print "A:\n", A
#print "a:\n", a
#print "b:\n", b
#print "c:\n", c
p_min = 0.2
mu_star = 1.5
zeta = 2.5
zeta_star = 1.5

k0 = 1e-3*np.matrix([-24.23, -0.2122, -10.5]).T
k1 = 1e-3*np.matrix([-40.13, -0.1396, -15.10]).T
k0_star = 1e-3*np.matrix([-11.96, -0.1111, -5.18]).T
k1_star = 1e-3*np.matrix([-83.35, -0.0477, -33.86]).T

##### Entwurf parameter #####
beta=2 # beta >=1 !

In [ ]:
###########################
# Hydraulischer Aktor     #
###########################

A0 = np.matrix([[0,   1,       0],
                [-10, -1.167, 25],
                [0,   0,    -0.8]])
print "Eigenvalues: {}".format(LA.eigvals(A0))
#a = -A[-1,:].T ### !!!!!
#print a
b0 = np.matrix([[0],[0],[2.4]])
c0 = np.matrix([1, 0, 0])
d0 = np.matrix([0])
u_max = 10.5
n = 3

X00 = [np.matrix([-20.0, -10.0, -10.0]).T,
       np.matrix([-20.0, -10.0, 10.0]).T,
       np.matrix([-20.0,  10.0, -10.0]).T,
       np.matrix([20.0,  -10.0, 10.0]).T,
       np.matrix([20.0,  -10.0, -10.0]).T,
       np.matrix([20.0,   10.0, 10.0]).T]

#print "A:\n", A
#print "a:\n", a
#print "b:\n", b
#print "c:\n", c

# Convert to Normalform
(A1, b1, c1, d1), T1, Q1 = optim_tools.get_Steuerungsnormalform(A0, b0, c0.T, d0)
a1 = -A1[-1][:].T #!!!!
print "T1:\n", T1

# Convert to Normalform
ss, T2 = con.canonical_form(con.ss(A0, b0, c0, d0), form='reachable')

assert np.allclose(T1*X00[1],
                   optim_tools.reverse_x_order(T2*X00[1])),\
"own Steuerungsnormalform Transformation not equal python control version"
#print "x_r1:\n", T1*X00[1]
#print "x_r2(backwards):\n", optim_tools.reverse_x_order(T2*X00[1])

A = optim_tools.reverse_x_order(np.matrix(ss.A))
a = -A[-1][:].T #!!!!

b = optim_tools.reverse_x_order(np.matrix(ss.B))
c = optim_tools.reverse_x_order(np.matrix(ss.C))
d = optim_tools.reverse_x_order(np.matrix(ss.D)) # == 0!

print "A:\n", A
assert np.allclose(A, A1)
print "a:\n", a
assert np.allclose(a, a1)
print "b:\n", b
assert np.allclose(b, b1)
print "c:\n", c
assert np.allclose(c, c1.T)

X0 = [T1.dot(x0) for x0 in X00]
print "X0:\n", X0
#print "A1:\n", A1
#print "a1:\n", a1
#print "b1:\n", b1
#print "c1:\n", c1