<s> """<STR_LIT>""" <EOL> from __future__ import division , print_function , absolute_import <EOL> import numpy as np <EOL> from numpy . linalg import inv , LinAlgError <EOL> from . basic import solve <EOL> from . lapack import get_lapack_funcs <EOL> from . decomp_schur import schur <EOL> from . special_matrices import kron <EOL> __all__ = [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , <EOL> '<STR_LIT>' , '<STR_LIT>' ] <EOL> def solve_sylvester ( a , b , q ) : <EOL> """<STR_LIT>""" <EOL> r , u = schur ( a , output = '<STR_LIT>' ) <EOL> s , v = schur ( b . conj ( ) . transpose ( ) , output = '<STR_LIT>' ) <EOL> f = np . dot ( np . dot ( u . conj ( ) . transpose ( ) , q ) , v ) <EOL> trsyl , = get_lapack_funcs ( ( '<STR_LIT>' , ) , ( r , s , f ) ) <EOL> if trsyl is None : <EOL> raise RuntimeError ( '<STR_LIT>' <EOL> '<STR_LIT>' ) <EOL> y , scale , info = trsyl ( r , s , f , tranb = '<STR_LIT:C>' ) <EOL> y = scale * y <EOL> if info < <NUM_LIT:0> : <EOL> raise LinAlgError ( "<STR_LIT>" <EOL> "<STR_LIT>" % ( - info , ) ) <EOL> return np . dot ( np . dot ( u , y ) , v . conj ( ) . transpose ( ) ) <EOL> def solve_lyapunov ( a , q ) : <EOL> """<STR_LIT>""" <EOL> return solve_sylvester ( a , a . conj ( ) . transpose ( ) , q ) <EOL> def _solve_discrete_lyapunov_direct ( a , q ) : <EOL> """<STR_LIT>""" <EOL> lhs = kron ( a , a . conj ( ) ) <EOL> lhs = np . eye ( lhs . shape [ <NUM_LIT:0> ] ) - lhs <EOL> x = solve ( lhs , q . flatten ( ) ) <EOL> return np . reshape ( x , q . shape ) <EOL> def _solve_discrete_lyapunov_bilinear ( a , q ) : <EOL> """<STR_LIT>""" <EOL> eye = np . eye ( a . shape [ <NUM_LIT:0> ] ) <EOL> aH = a . conj ( ) . transpose ( ) <EOL> aHI_inv = inv ( aH + eye ) <EOL> b = np . dot ( aH - eye , aHI_inv ) <EOL> c = <NUM_LIT:2> * np . dot ( np . dot ( inv ( a + eye ) , q ) , aHI_inv ) <EOL> return solve_lyapunov ( b . conj ( ) . transpose ( ) , - c ) <EOL> def solve_discrete_lyapunov ( a , q , method = None ) : <EOL> """<STR_LIT>""" <EOL> a = np . asarray ( a ) <EOL> q = np . asarray ( q ) <EOL> if method is None : <EOL> if a . shape [ <NUM_LIT:0> ] >= <NUM_LIT:10> : <EOL> method = '<STR_LIT>' <EOL> else : <EOL> method = '<STR_LIT>' <EOL> meth = method . lower ( ) <EOL> if meth == '<STR_LIT>' : <EOL> x = _solve_discrete_lyapunov_direct ( a , q ) <EOL> elif meth == '<STR_LIT>' : <EOL> x = _solve_discrete_lyapunov_bilinear ( a , q ) <EOL> else : <EOL> raise ValueError ( '<STR_LIT>' % method ) <EOL> return x <EOL> def solve_continuous_are ( a , b , q , r ) : <EOL> """<STR_LIT>""" <EOL> try : <EOL> g = inv ( r ) <EOL> except LinAlgError : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' ) <EOL> g = np . dot ( np . dot ( b , g ) , b . conj ( ) . transpose ( ) ) <EOL> z11 = a <EOL> z12 = - <NUM_LIT:1.0> * g <EOL> z21 = - <NUM_LIT:1.0> * q <EOL> z22 = - <NUM_LIT:1.0> * a . conj ( ) . transpose ( ) <EOL> z = np . vstack ( ( np . hstack ( ( z11 , z12 ) ) , np . hstack ( ( z21 , z22 ) ) ) ) <EOL> s , u , _ = schur ( z , sort = '<STR_LIT>' ) <EOL> ( m , n ) = u . shape <EOL> u11 = u [ <NUM_LIT:0> : m // <NUM_LIT:2> , <NUM_LIT:0> : n // <NUM_LIT:2> ] <EOL> u21 = u [ m // <NUM_LIT:2> : m , <NUM_LIT:0> : n // <NUM_LIT:2> ] <EOL> u11i = inv ( u11 ) <EOL> return np . dot ( u21 , u11i ) <EOL> def solve_discrete_are ( a , b , q , r ) : <EOL> """<STR_LIT>""" <EOL> try : <EOL> g = inv ( r ) <EOL> except LinAlgError : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' ) <EOL> g = np . dot ( np . dot ( b , g ) , b . conj ( ) . transpose ( ) ) <EOL> try : <EOL> ait = inv ( a ) . conj ( ) . transpose ( ) <EOL> except LinAlgError : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' ) <EOL> z11 = a + np . dot ( np . dot ( g , ait ) , q ) <EOL> z12 = - <NUM_LIT:1.0> * np . dot ( g , ait ) <EOL> z21 = - <NUM_LIT:1.0> * np . dot ( ait , q ) <EOL> z22 = ait <EOL> z = np . vstack ( ( np . hstack ( ( z11 , z12 ) ) , np . hstack ( ( z21 , z22 ) ) ) ) <EOL> s , u , _ = schur ( z , sort = '<STR_LIT>' ) <EOL> ( m , n ) = u . shape <EOL> u11 = u [ <NUM_LIT:0> : m // <NUM_LIT:2> , <NUM_LIT:0> : n // <NUM_LIT:2> ] <EOL> u21 = u [ m // <NUM_LIT:2> : m , <NUM_LIT:0> : n // <NUM_LIT:2> ] <EOL> u11i = inv ( u11 ) <EOL> return np . dot ( u21 , u11i ) </s>
