In [ ]:
import numpy as np
import control as con
from numpy import linalg as LA
import cvxpy
import optim_tools #own file with helper
import boris_tools
import matplotlib.pyplot as plt
%matplotlib inline
In [ ]:
print cvxpy.installed_solvers()
In [ ]:
##############################
# Boris Diss Reactor #
##############################
A = np.matrix([[1, 0],
[0, -0.5]])
b = np.matrix([[1],[-0.5]])
c = np.matrix([1, 0])
d = np.matrix([0])
u_max = 1
n = 2
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 [ ]:
### Design Sättigungsregler mittels konvexer Hülle (A.3)
# Variables (Convex)
# Name in A.3 | Name in Program
# Q = P^⁻1 | = Q
# z = Ql | = z0
# z* = Ql* | = z1
# Variables (Quasi convex)
# gamma | = g
# Parameter
# mu | = mu (Designparameter)
# X0 | = X0 = [x0,...]
# A | = A
# b | = b
# Initialize
n = len(b) # get dim of system
# Define Variables
Q = cvxpy.Semidef(n, n) #symmetric and positive semidefinite
#Q = cvxpy.Variable(n, n) # Semidef could go as an additional constraint as well, thereby no need fo Q to be symmetric
#const_sdQ = Q >> 0 # semidef but not (necessarily) symmetric
z0 = cvxpy.Variable(n)
z1 = cvxpy.Variable(n)
# Bisection parameter
g = cvxpy.Parameter(sign='positive')
#g.value = g_val # TODO set by bisection loop (function?)
mu = cvxpy.Parameter(sign='positive')
#m.value = 1 # TODO: mu* >=1
# Define Constraints
# (A.10)
const_A10 = [cvxpy.bmat([[Q, X0[i]],
[X0[i].T, 1 ]]) >> 0
for i in range(0, len(X0))]
# (A.11)
const_A11 = cvxpy.bmat([[Q, z0],
[z0.T, 1 ]]) >> 0
# (A.12)
const_A12 = cvxpy.bmat([[Q, z1 ],
[z1.T, mu**2]]) >> 0
# (A.13)
const_A13 = Q*A + A.T*Q - b*z0.T - z0*b.T << 0 # This constraint is strict definit
# (A.14)
const_A14 = Q*A + A.T*Q - b*z1.T - z1*b.T << -2*g*Q # This constraint is strict definit
#const_A14 = Q*A + A.T*Q - b*z1.T - z1*b.T << cvxpy.bmat([[-2*g, 0, 0],
# [0, -2*g, 0],
# [0, 0, -2*g]])*Q # This constraint is strict definit
# Collect all constraints
constraints_A15 = []
constraints_A15.extend(const_A10) ##!! Beware of the "extend" if input is array
constraints_A15.append(const_A11)
constraints_A15.append(const_A12)
constraints_A15.append(const_A13)
constraints_A15.append(const_A14)
# Feasibility for bisection:
obj_A15 = cvxpy.Minimize(0)
# Stronger requirements for bisection? -> probably resulting in higher iterations, thus better results: TODO: Fix bisect!
# obj_alt = cvxpy.Maximize(cvxpy.log_det(Q)) # Identical to geo_mean (in term of convexity and result)
# Form and solve problem.
prob_A15 = cvxpy.Problem(obj_A15, constraints_A15)
In [ ]:
%%time
# Test bisection
mu.value = 1 # TODO: mu* >=1
[[Q_o, z0_o, z1_o], g_o] = optim_tools.bisect_max(0, None, prob_A15, g, [Q, z0, z1], bisect_verbose=False,
solver=cvxpy.CVXOPT)
print "Sättigungsregler mittels konvexer Hülle -> Bisection(Max. Abklingrate)"
print "Bisection Param:", g_o
# Output
print "Status:", prob_A15.status
Q_A15 = Q_o
P_A15 = LA.inv(Q_o)
l0_A15 = P_A15.dot(z0_o)
l1_A15 = P_A15.dot(z1_o)
print "Q:\n", Q_A15
print "P:\n", P_A15
print "z:\n", z0_o
print "z*:\n", z1_o
print "l:\n", l0_A15
print "l*:\n", l1_A15
print
In [ ]:
# Werte nach Tabelle 4.5
zeta0 = [1.1, 1.5, 2.5, 4, 7]
zeta1 = [1.5, 2, 3, 4, 5, 8]
mu1 = [1.5, 3, 5]
pmin = [1.0/10, 1.0/20]
In [ ]:
# Selection according to underlined results Tab 4.5
#TODO: Loop over four arrays
i_zeta0 = 2
i_zeta1 = 4
i_mu1 = 0
i_pmin = 1
print pmin[i_pmin]
print mu1[i_mu1]
print zeta0[i_zeta0]
print zeta1[i_zeta1]
In [ ]:
print "Entering step 1"
In [ ]:
%%time
# Perform bisection on g for mu=1
mu.value = 1
[[Q_o, z0_o, z1_o], g_o] = optim_tools.bisect_max(0, None, prob_A15, g, [Q, z0, z1], bisect_verbose=False,
solver=cvxpy.CVXOPT)
print "-----------RESULTS-----------"
print "Bisection Param:", g_o
# Output
print "Status:", prob_A15.status
# l*(mu=1 -> m0) -> l1m0
l1m0 = LA.inv(Q_o) * z1_o
print "l*(mu=1) =\n", l1m0
roots_m0_p1 = LA.eigvals(A-b*l1m0.T)
print "lambda^{hat}(p=1) = ", roots_m0_p1
print
In [ ]:
%%time
# Perform bisection on g for mu*>=1 (m1)
mu.value = mu1[i_mu1] # TODO: mu* comes from optimisation loop
[[Q_o, z0_o, z1_o], g_o] = optim_tools.bisect_max(0, None, prob_A15, g, [Q, z0, z1], bisect_verbose=False,
solver=cvxpy.CVXOPT)
print "-----------RESULTS-----------"
print "Bisection Param:", g_o
# Output
print "Status:", prob_A15.status
# l*(mu=1 -> m0) -> l1m0
l1m1 = LA.inv(Q_o) * z1_o
print "l*(mu=1) =\n", l1m1
roots_m1_p1 = LA.eigvals(A-b*l1m1.T)
print "lambda^{hat}(p=1) = ", roots_m1_p1
print
In [ ]:
roots_m0_pmin = zeta0[i_zeta0] * roots_m0_p1.real + 1j*roots_m0_p1.imag
#roots_m0_pmin = zeta0 * roots_m0_p1 # Shifting on both real and imaginary axis in zeta
print "lambda^{hat}(p=pmin) = ", roots_m0_pmin
In [ ]:
roots_m1_pmin = zeta1[i_zeta1] * roots_m1_p1.real + 1j*roots_m1_p1.imag
#roots_m0_pmin = zeta0 * roots_m0_p1 # Shifting on both real and imaginary axis in zeta
print "lambda^{hat}(p=pmin) = ", roots_m1_pmin
In [ ]:
# Get Polynomcoefficiants
# k0_i -> k_i
k0_0, k0_1 = boris_tools.k_explizit_Ab2(roots_m0_p1, roots_m0_pmin, pmin[i_pmin], A, b)
In [ ]:
# Plot poles from A^{hat}(p)
plt.axis([-3, 0, -1.5, 1.5])
boris_tools.plot_moving_poles(A, b, c, d, k0_0, k0_1, pmin=pmin[i_pmin])
In [ ]:
# Get Polynomcoefficiants
# k1_i -> k_{*,i}
k1_0, k1_1 = boris_tools.k_explizit_Ab2(roots_m1_p1, roots_m1_pmin, pmin[i_pmin], A, b)
In [ ]:
plt.axis([-3, 0, -1.5, 1.5])
boris_tools.plot_moving_poles(A, b, c, d, k1_0, k1_1, pmin[i_pmin])
print pmin[i_pmin]
In [ ]:
print "Entering step 2"
In [ ]:
# Transformations A.4
from scipy.special import comb as nchoosek # n Choose k (n ueber k)
# n muss be dim(Q) = dim(R1) = dim(A)
# manual approved with equations
def get_S_list(u_max, Q, z, a, n):
S_list = [cvxpy.bmat([[u_max**2 , z.T],
[z , Q]])] # S0 is different !
for i in range(n+1)[1:]:
#print i
#print -a[n-i]
q = Q[:, n-i] # Slicing, (n-i)th column of Q!
#print q
S_list.append(-a[n-i] * cvxpy.bmat([[0, q.T],
[q, np.zeros((n,n))]]))
return S_list
# As seen in A.4 jasniwiecz
# Intervalltransformation (Gleichung (4.36))
# p => [p_l, p_u] (p => [p_min, 1]) wird zu p^ => [-1,1] return neue Matrizen(!) mit diesen p^
#
# S_in ist Liste aus Faktormatrizen für Polynom in p
# S_out ist Liste aus Faktormatrizen für Polynom in p^
# a = (p_u - p_l)
# b = (p_u + p_l)
# manual approved with equations
def trans_S_list(S_in, pl=0.1, pu=1):
a = pu - pl
b = pu + pl
n = len(S_in) # Anzahl der Matrizen in S_in
S_out = [np.zeros(S_in[0].size)] * n
for j in range(0, n): # Bearbeite jede Matrix
for i in range(j, n): # Jeweils mit Summe der anderen
S_out[j] = S_out[j] + 2**(-i) * b**(i-j) * nchoosek(i, j) * S_in[i]
S_out[j] = a**j*S_out[j]
return S_out
def calc_S_Sum(S_list):
l = len(S_list) # number of matrizen in S_list
m = l-1 # Index of Q_m
n = S_list[0].size[0] # shape of each matrix, first element
if m is 0:
S_sum = cvxpy.bmat([[2*S_list[0], np.zeros(n)],
[np.zeros(n), np.zeros(n)]])
elif m is 1:
S_sum = cvxpy.bmat([[2*S_list[0], S_list[1]],
[S_list[1], np.zeros(n)]])
else: # e.g. m is 2 or more
S_sum = cvxpy.bmat([[2*S_list[0], S_list[1]],
[S_list[1], 2*S_list[2]]])
for i1 in range(3, l, 2):
S_new_col = cvxpy.vstack(np.zeros((((i1+1)/2-1)*n, n)), S_list[i1])
if i1 is m:
S_new_row = cvxpy.hstack(np.zeros((n, ((i1+1)/2-1)*n)), S_list[i1], np.zeros((n,n)))
else:
S_new_row = cvxpy.hstack(np.zeros((n, ((i1+1)/2-1)*n)), S_list[i1], 2*S_list[i1+1])
S_sum = cvxpy.bmat([[S_sum, S_new_col],
[S_new_row]])
# Beware, there is NO minus compared with Dilyanas version
S_sum = 0.5*S_sum
return S_sum
# CJ = "Selection matrizes" of (31), l=dimension of P and G
def calc_lmi_cond(S_sum, n):
k = S_sum.size[1] / n
J = np.hstack([np.zeros((n*(k-1), n)), np.eye(n*(k-1))])
C = np.hstack([np.eye(n*(k-1)), np.zeros((n*(k-1), n))])
CJ = np.vstack([C, J])
l = n*(k-1)
return CJ, l
In [ ]:
###############################
# Optimisation problem (4.40) #
###############################
# Define Variables
R0 = cvxpy.Semidef(n, n)
R1 = cvxpy.Semidef(n, n)
G0 = cvxpy.Variable(n, n) # G -> skew
G1 = cvxpy.Variable(n, n) # G_* -> skew
D0 = cvxpy.Variable(n, n) # D -> sym
D1 = cvxpy.Variable(n, n) # D_* -> sym
# Define Parameters
pm = cvxpy.Parameter(sign='positive') # pmin
# TODO: pmin from optimazation loop
pm.value = pmin[i_pmin]
alpha = 1-pm
beta = 1+pm
# TODO: Check if Â_0 and Â_1 are correct (find proof)
A0_0 = A-b*k0_0.T # A^{hat}_0
A0_1 = b*k0_1.T # A^{hat}_0
#A0_1 = A-b*k0_1.T # A^{hat}_1
T0_0 = A0_1.T*R1 + R1*A0_1 # Theta_0
T0_1 = A0_0.T*R1 + R1*A0_0 + A0_1.T*R0 + R0*A0_1 # Theta_1
T0_2 = A0_0.T*R0 + R0*A0_0 # Theta_2
# TODO: Check if Â_{*,0} and Â_{*,1} are correct (find proof)
A1_0 = A-b*k1_0.T # A^{hat}_{*, 0}
A1_1 = b*k1_1.T # A^{hat}_{*, 0}
#A1_1 = A-b*k1_1.T # A^{hat}_{*, 1}
T1_0 = A1_1.T*R1 + R1*A1_1 # Theta_{*,0}
T1_1 = A1_0.T*R1 + R1*A1_0 + A1_1.T*R0 + R0*A1_1 # Theta_{*,1}
T1_2 = A1_0.T*R0 + R0*A1_0 # Theta_{*,2}
# Define constraints
# Constraints on helper variables
constraint_D0 = D0 == D0.T # symmetry
constraint_D1 = D1 == D1.T # symmetry
constraint_G0 = G0 == -G0.T # skew symmetry
constraint_G1 = G1 == -G1.T # skew symmetry
# Constraints (4.39)
constraint_439a = R0 >> 0 # Feasible by definition (R0 is semidef)
constraint_439b = R0 + R1 >> 0 # Feasible by definition (R1 is semidef)
constraint_439c = cvxpy.bmat([[-D0 , G0],
[G0.T, D0]]) - \
cvxpy.bmat([[T0_0 + 0.5*beta*T0_1 + 0.25*beta**2*T0_2, 0.25*alpha*(T0_1 + beta*T0_2)],
[0.25*alpha*(T0_1 + beta*T0_2), 0.25*alpha**2*T0_2 ]]) >> 0
constraints_439d = [1.0 - X0[i].T*(R0 + R1)*X0[i] > 0 for i in range(0, len(X0))]
constraint_439e = cvxpy.bmat([[1, k0_0.T],
[k0_0, R0]]) + \
1.0/pm * cvxpy.bmat([[0, k0_1.T],
[k0_1, R1]]) >> 0
constraint_439e_ = cvxpy.bmat([[1, k0_0.T+1.0/pm*k0_1.T],
[k0_0+1.0/pm*k0_1, R0+1.0/pm*R1]]) >> 0
constraint_439f = cvxpy.bmat([[1, k0_0.T],
[k0_0, R0]]) + \
cvxpy.bmat([[0, k0_1.T],
[k0_1, R1]]) >> 0
constraint_439f_ = cvxpy.bmat([[1, k0_0.T+k0_1.T],
[k0_0+k0_1, R0+R1]]) >> 0
constraint_439g = cvxpy.bmat([[-D1 , G1],
[G1.T, D1]]) - \
cvxpy.bmat([[T1_0 + 0.5*beta*T1_1 + 0.25*beta**2*T1_2, 0.25*alpha*(T1_1 + beta*T1_2)],
[0.25*alpha*(T1_1 + beta*T1_2), 0.25*alpha**2*T1_2 ]]) >> 0
constraints_439 = [constraint_D0, # sym
constraint_D1, # sym
constraint_G0, # skew
constraint_G1, # skew
constraint_439a,
constraint_439b]
constraints_439.append(constraint_439c)
constraints_439.extend(constraints_439d) ##!! Beware of the "extend" if input is array
constraints_439.append(constraint_439e)
constraints_439.append(constraint_439f)
constraints_439.append(constraint_439g)
# Form objective.
obj_440 = cvxpy.Minimize(cvxpy.trace(R0+1.0/pm*R1))
# Form and solve problem.
prob_440 = cvxpy.Problem(obj_440, constraints_439)
# TODO: Make complete or remove
def check_cons():
try:
print "R0", R0.value
LA.cholesky(R0.value)
print "R1", R1.value
LA.cholesky(R1.value + R0.value)
print "\n\n-----------------------------------\nSolution ok!"
except:
print "\n\n-----------------------------------\nSolution NOT ok!"
pass
try:
prob_440.solve(solver=cvxpy.SCS, verbose=False, max_iters=5000) # Returns the optimal value.
print "SCS status:", prob_440.status
check_cons()
except:
print "SCS failed"
try:
prob_440.solve(solver=cvxpy.CVXOPT, verbose=False) # Returns the optimal value.
print "CVXOPT status:", prob_440.status
check_cons()
except:
print "CVXOPT failed"
try:
prob_440.solve(solver=cvxpy.MOSEK, verbose=False) # Returns the optimal value.
print "MOSEK status:", prob_440.status
check_cons()
except:
print "MOSEK failed"
#print "optimal value", prob_440.value
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: