QR decomposition

schur_qr.py

import copy
from casadi import *
import numpy as np

def qr_eigen(A, iterations=100):
    pQ = SX.eye(A.size1())
    X = copy.deepcopy(A)
    for _ in range(iterations):
        Q, R = qr(X)  # QR decomposition in CasADi
        pQ = pQ @ Q # Update eigenvector matrix for next iteration
        X = R @ Q  # Update eigenvalue matrix for next iteration
    return diag(X), pQ


n_X = 50
N = 100
KSx = SX.sym('KSx', n_X,n_X)
# KMx = MX.sym('KSx', 10,10)

w_sym , v_sym = qr_eigen(KSx, 30)

eigFun = Function('WV', [KSx], [w_sym])

# eigMX = eigFun(KMx)

opti =Opti()

K = opti.variable(n_X,n_X)
Dx = opti.parameter(n_X, N)
Dy = opti.parameter(n_X, N)

opti.minimize(sumsqr(K@Dx - Dy))

K_eig = eigFun(K)

for k in range(K_eig.size()[0]):
    opti.subject_to(opti.bounded(0.0, norm_2(K_eig), 1.0))

# SOLVER SETTINGS
# Setting without hessian
# opti.solver('ipopt', {}, {'hessian_approximation': 'limited-memory', "max_iter":6000, "limited_memory_update_type": "sr1"})
opti.solver('sqpmethod',{"qpsol":"qrqp","convexify_strategy": "regularize", "hessian_approximation": "GN", "print_iteration":True,"print_time":True,"print_status":True,"print_header":False,"max_iter":10000,"qpsol_options": {"print_iter":False,"print_header":False}})

# Setting with hessian
# opti.solver('sqpmethod',{"qpsol":"qrqp","convexify_strategy": "regularize", "print_iteration":True,"print_time":True,"print_status":True,"print_header":False,"max_iter":10000,"qpsol_options": {"print_iter":False,"print_header":False}})
# opti.solver("ipopt", {})


J = opti.to_function('K_J',[Dx, Dy], [K])