In [ ]:
# intervalltransformation
# i = i2
# j = i1
import numpy as np
from scipy.special import comb as nchoosek # n Choose k (n ueber k)
def transQ(Q_in, a,b):
#check array
if not isinstance(Q_in, list):
print "list input expected"
# m,n
m = len(Q_in)
n = (Q_in[0]).shape[0]
Q_out = Q_in
for i1 in range (0,m):
#Q_out[i1,:] = Q_in[i1]
Q_out[i1] = Q_in[i1]*2**(-i1)
for i2 in range (i1+1, m):
#print Q_out[i1]
Q_out[i1] = Q_out[i1] + (2**(-i2))*(b**(i2-i1))*nchoosek(i2,i1)*Q_in[i2];
Q_out[i1] = (a**i1)*Q_out[i1]
# nested loop
return Q_out
In [ ]:
transQ([np.array([[1, 1],[1, 1]]),
np.array([[2, 2],[2, 2]]),
np.array([[3, 3],[3, 3]])], 1, 1)
In [ ]:
# intervalltransformation
import numpy as np
def calcQsum(Q_in):
if not isinstance(Q_in, list):
print "list input expected"
m = len(Q_in)
print m
n = (Q_in[0]).shape[0]
if m == 1:
#Q_sum = np.matrix([[2*Q_in[0], np.zeros((n,n))],
# [np.zeros((n,n)), np.zeros((n,n))]])
#np.column_stack((a,b))
Q_sum = np.matrix(np.row_stack((np.column_stack((2*Q_in[0], np.zeros((n,n)))),
np.column_stack((np.zeros((n,n)), np.zeros((n,n)))))))
elif m == 2:
#Q_sum = np.matrix([[2*Q_in[0], Q_in[1,:]],
# [Q_in[1,:], np.zeros((n,n))]])
Q_sum = np.matrix(np.row_stack((np.column_stack((2*Q_in[0], Q_in[1])),
np.column_stack((Q_in[1], np.zeros((n,n)))))))
else:
#print "wthin"
#Q_sum = np.matrix([[2*Q_in[0], Q_in[1,:]],
# [Q_in[1,:], 2*Q_in[2,:]]])
Q_sum = np.matrix(np.row_stack((np.column_stack((2*Q_in[0], Q_in[1])),
np.column_stack((Q_in[1], 2*Q_in[2])))))
#print "wth"
#print Q_sum
loopRange = np.arange(3,m,2)
for i1 in loopRange:
if i1 != (m-1):
#print Q_sum
#print Q_in[i1+1]
Q_sum = np.matrix(np.row_stack((np.column_stack((Q_sum,
np.row_stack(([np.zeros((int((i1+1)*0.5-1)*n,n)),
Q_in[i1]])))),
np.column_stack((np.zeros((n,int((i1+1)*0.5-1)*n)),
Q_in[i1],
2*Q_in[i1+1])))))
else:
print i1*0.5-1
#print n
#print np.zeros((int((i1+1)*0.5-1)*n,n))
#print Q_in[i1]
Q_sum = np.matrix(np.row_stack((np.column_stack((Q_sum,
np.row_stack((np.zeros((int((i1+1)*0.5-1)*n,n)), Q_in[i1])))),
np.column_stack((np.zeros((n,(int((i1+1)*0.5-1)*n))),
Q_in[i1],
np.zeros((n,n)))))))
Q_sum = (0.5)*Q_sum;
return Q_sum
In [ ]:
testarr = [np.array([[1, 1],
[1, 1]]),
np.array([[2, 2],
[2, 2]]),
np.array([[3, 3],
[3, 3]]),
np.array([[4, 4],
[4, 4]]),
np.array([[5, 5],
[5, 5]])]
calcQsum(testarr[0:5])
In [ ]:
# Selection stratergy
# plotting Boundary with given p
# g(p,x) = e(p).x^T.R.x-1
# p: (0 ,1]
# R:
# a_head: coefficiences os fixed-"zulegenden" characteristic
# Polynomials des closed control-loop with p = 1
# a: coefficiences of characteristic "Strecke"
import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt
from sympy.solvers import solve
from sympy import MatrixSymbol, Matrix
def boundaryPlotting_LyapunovFunction(p, R, a_head, a, stepJump):
n = np.size(a_head,0)
p_array = np.zeros(n)
for i in range(0, n) :
p_array[i] = np.power(p, n-i)
D_p = np.diag(p_array)
K_p = inv(D_p).dot(a_head) - a
R_p = inv(D_p).dot(R).dot(inv(D_p))
e_p = np.transpose(K_p).dot(inv(R_p)).dot(K_p)
# print "R1 = \n",R
# print "a_head = \n",a_head
# print "a = \n",a
# print "n = \n",n
# print "D_p = \n",D_p
# print "K_p = \n",K_p
# print "R_p = \n",R_p
# print "e_p = \n",e_p
result = []
ind = 0
# for x1 in np.arange(-0.5,0.5,0.004):
for x1 in np.arange(-0.5,0.5,stepJump):
for x2 in np.arange(-0.5,0.5,stepJump):
X = np.matrix([[x1],[x2]])
#print X.T
g_px = np.asscalar(e_p)*(X.T).dot(R_p).dot(X)-1
# print g_px
if (abs(g_px) <= 0.0002):
# plt.plot(x1, x2, 'ro',label='closed')
result.append(([x1,x2]))
ind = ind +1
# plt.show()
return result
In [ ]:
import matplotlib.patches as mpatches
R1 = np.matrix([[4, 1],
[1, 2]])
a_head = np.matrix([[2],
[2]])
a = np.matrix([[0],
[0]])
p1 = boundaryPlotting_LyapunovFunction(0.99999, R1, a_head, a, 0.002)
p2 = boundaryPlotting_LyapunovFunction(0.9, R1, a_head, a, 0.002)
p3 = boundaryPlotting_LyapunovFunction(0.6, R1, a_head, a, 0.002)
p4 = boundaryPlotting_LyapunovFunction(0.4, R1, a_head, a, 0.002)
In [ ]:
def arrangeListToPlot(p1):
#print len(p1)
r = max(p1, key=lambda item: item[0]) # most left of eclipse
l = min(p1, key=lambda item: item[0]) # most right of eclipse
u = max(p1, key=lambda item: item[1]) # most up of eclipse
d = min(p1, key=lambda item: item[1]) # most down of eclipse
print r
print l
print u
print d
meanp1 = np.mean(p1,0)
print meanp1
p1u = []
p1d = []
p1ul = []
p1ur = []
p1dl = []
p1dr = []
# seperate to bigger and smaller or equal mean of second parameter
#for item in p1:
# if (item[1] > meanp1[1]):
# p1u.append(item)
# else:
# p1d.append(item)
# seperate to bigger and smaller or equal mean of first parameter
#for item in p1u:
# if (item[0] > meanp1[0]):
# p1ur.append(item)
# else:
# p1ul.append(item)
# seperate to bigger and smaller or equal mean of first parameter
#for item in p1d:
# if (item[0] > meanp1[0]):
# p1dr.append(item)
# else:
# p1dl.append(item)
for item in p1:
if (item[1] >= l[1] and item[0] <= u[0]):
p1ul.append(item)
elif (item[1] >= r[1] and item[0] >= u[0]):
p1ur.append(item)
elif (item[1] <= r[1] and item[0] >= d[0]):
p1dr.append(item)
elif (item[1] <= l[1] and item[0] <= d[0]):
p1dl.append(item)
#print len(p1ul)
#print len(p1ur)
#print len(p1dl)
#print len(p1dr)
list.sort(p1ul,key=lambda l:l[1])
list.sort(p1ur,key=lambda l:l[1], reverse=True)
list.sort(p1dr,key=lambda l:l[1], reverse=True)
list.sort(p1dl,key=lambda l:l[1])
tmp = []
tmp.extend(p1ul)
tmp.extend(p1ur)
tmp.extend(p1dr)
tmp.extend(p1dl)
tmp.extend(p1ul)
#print p1dr
print len(tmp)
return tmp
In [ ]:
tmp1 = arrangeListToPlot(p1)
tmp2 = arrangeListToPlot(p2)
tmp3 = arrangeListToPlot(p3)
tmp4 = arrangeListToPlot(p4)
p1x = (np.asarray(tmp1))[:,0]
p1y = (np.asarray(tmp1))[:,1]
p2x = (np.asarray(tmp2))[:,0]
p2y = (np.asarray(tmp2))[:,1]
p3x = (np.asarray(tmp3))[:,0]
p3y = (np.asarray(tmp3))[:,1]
p4x = (np.asarray(tmp4))[:,0]
p4y = (np.asarray(tmp4))[:,1]
#[x for x in X if P(x)]
#p1u = [item for item in p1 if item]
plt.plot(p1x,p1y, '-o',label='closed')
plt.plot(p2x,p2y, '-o',label='closed')
plt.plot(p3x,p3y, '-o',label='closed')
plt.plot(p4x,p4y, '-o',label='closed')
plt.show()
#tmp1
LORENZ RETRY
In [ ]:
import numpy as np
from numpy import linalg as LA
# Harmonischer Oszillator
A = np.matrix([[0, 1],
[-1, 0]])
b = np.matrix([0, 1]).T
R1 = np.matrix([[1, 0.975],
[0.9750, 2.389]])
X0 = [np.matrix([1.5, 1.5]),
np.matrix([1.5, -1.5]),
np.matrix([-1.5, 1.5]),
np.matrix([-1.5, -1.5])]
a = np.matrix([1, 0]).T
a_hat = np.matrix([2, np.sqrt(6)]).T
In [ ]:
def _D_inv(v, n):
return np.diag([v**-x for x in range(n, 0, -1)])
# k(v) (4.5)
def get_k(v, a, a_hat):
try:
v = v.squeeze() # This is weird, needed for minimize for some reason
except:
pass
D_inv = _D_inv(v, len(a))
#print D_inv
k = D_inv.dot(a_hat) - a
return k
#print "Test k\n", get_k(0.2, np.matrix([124, 48, 12]).T, np.matrix([1, 2, 3]).T)
#assert(np.allclose(-a+a_hat, get_k(1, a, a_hat)))
# Fixing a and a_hat for convinience
# func_k = lambda v: get_k(v, a, a_hat)
# G(v) (4.4)
def get_g(v, x, R1, u_max, a, a_hat):
try:
v = v.squeeze() # This is weird, needed for minimize for some reason
except:
pass
D_inv = _D_inv(v, len(x))
R = D_inv.dot(R1).dot(D_inv) # R(v) = D^⁻1 * R1 * D^-1
k = get_k(v, a, a_hat) # k(v)
e = (u_max**(-2)) * (k.T.dot(LA.inv(R)).dot(k))
#assert e < 1.0
g = e*(x.T.dot(R).dot(x)) - 1.0
#assert g <= 0, "g = {} > 0".format(g)
# Update 2016: As of python 3.5, there is a new matrix_multiply symbol, @:
# g = x' @ D^-1 @ R1 @ D^-1 @ x - 1.0
return g
In [ ]:
from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt
#points = np.random.rand(30, 2) # 30 random points in 2-D
#hull = ConvexHull(points)
#print X0
#print len(X0)
# Convert boundaries
poi = X0[0]
for i in range(1, len(X0)):
poi = np.vstack((poi, X0[i]))
#print poi
# get convex hull
hull = ConvexHull(poi)
#plot it
plt.plot(poi[:,0], poi[:,1], 'o')
for simplex in hull.simplices:
plt.plot(poi[simplex, 0], poi[simplex, 1], 'k-')
# boundaries ready
# get elipsoid
# g = e*x.T*R*x -1
area = np.arange(-10, 10, .0001)
v = 1
e_poi = np.array([0 , 0])
for x1, x2 in zip(area, area):
if np.absolute(get_g(v, np.matrix([x1, x2]).T, R1, 10, a, a_hat)-1)<0.1:
e_poi = np.vstack((e_poi, [x1, x2]))
#print e_poi
# get convex hull
hull2 = ConvexHull(e_poi)
#plot it
plt.plot(e_poi[:,0], e_poi[:,1], 'o')
for simplex in hull2.simplices:
plt.plot(e_poi[simplex, 0], e_poi[simplex, 1], 'r-')
plt.show()
In [ ]:
plt.plot(e_poi[:,0], e_poi[:,1], 'o')
plt.show()
In [ ]:
In [ ]:
#plt.plot(points[hull.vertices,0], points[hull.vertices,1], 'r--', lw=2)
#plt.plot(points[hull.vertices[0],0], points[hull.vertices[0],1], 'ro')
In [ ]: