Code Reference

Projected Riccati equations

proj_ric_utils.comp_proj_lyap_res_norm(Z, amat=None, mmat=None, wmat=None, jmat=None, umat=None, vmat=None, Sinv=None)

compute the squared f norm of projected lyap residual

res = Pt*[ Ft*ZZt*M + Mt*ZZt*M + W*Wt ]*P

proj_ric_utils.compress_ZQR(Z, kmax=None, shplot=False)

routine that compresses the columns Z by means of rank revealing QR

such that it ZZ.T is still well approximated

proj_ric_utils.compress_Zsvd(Z, k=None, thresh=None, shplot=False, verbose=True)

routine that compresses the columns Z by means of a truncated SVD

such that it ZZ.T is still well approximated

proj_ric_utils.get_mTzzTtb(MT, Z, tB, output=None)

compute the left factor of the lyapunov coefficient

related to the linearization

proj_ric_utils.proj_alg_ric_newtonadi(mmat=None, amat=None, jmat=None, bmat=None, wmat=None, z0=None, mtxoldb=None, transposed=False, nwtn_adi_dict={'nwtn_max_steps': 14, 'nwtn_upd_reltol': 1e-08, 'adi_max_steps': 150, 'adi_newZ_reltol': 1e-05}, **kw)

solve the projected algebraic ricc via newton adi

M.T*X*A + A.T*X*M - M.T*X*B*B.T*X*M + J(Y) = -WW.T

JXM = 0 and M.TXJ.T = 0

If mtxb is given, (e.g. as the feedback computed in a previous step of a Newton iteration), the coefficient matrix with feedback

A.T <- A.T - mtxb*b

is considered

proj_ric_utils.pymess_dae2_newtonadi(mmat=None, amat=None, jmat=None, bmat=None, wmat=None, z0=None, mtxoldb=None, transposed=False, nwtn_adi_dict={'nwtn_max_steps': 14, 'nwtn_upd_reltol': 1e-08, 'adi_max_steps': 150, 'adi_newZ_reltol': 1e-05}, **kw)

solve the projected algebraic ricc via newton adi

M.T*X*A + A.T*X*M - M.T*X*B*B.T*X*M + J(Y) = -WW.T

JXM = 0 and M.TXJ.T = 0

If mtxb is given, (e.g. as the feedback computed in a previous step of a Newton iteration), the coefficient matrix with feedback

A.T <- A.T - mtxb*b

is considered

proj_ric_utils.solve_proj_lyap_stein(amat=None, jmat=None, wmat=None, mmat=None, umat=None, vmat=None, transposed=False, adi_dict={'adi_max_steps': 150, 'adi_newZ_reltol': 1e-08}, nwtn_adi_dict=None, **kw)

approximates the solution X to the projected lyap equation

[A-UV].T*X*M + M.T*X*[A-UV] + J.T*Y*M + M.T*Y.T*J = -W*W.T

J*X*M = 0 and M.T*X*J.T = 0

by considering the equivalent Stein eqns and computing the first members of the series converging to X

We use the SMW formula: (A-UV).-1 = A.-1 + A.-1*U [ I - V A.-1 U].-1 A.-1

for the transpose:

(A-UV).-T = A.-T + A.-T*Vt [ I - Ut A.-T Vt ].-1 A.-T

= (A.T - Vt Ut).-1

see numOptAff.pdf

proj_ric_utils.solve_stst_feedbacknthrough(amat=None, mmat=None, jmat=None, bmat=None, cmat=None, fv=None, fl=None, fg=None, fl2=None, nwtn_adi_dict={'nwtn_max_steps': 14, 'nwtn_upd_reltol': 1e-08, 'adi_max_steps': 150, 'adi_newZ_reltol': 1e-05})

solve for the stabilizing feedback gain and the feedthrough

for the linear time invariant case

Linear algebra utilities

lin_alg_utils.app_prj_via_sadpnt(amat=None, jmat=None, rhsv=None, jmatT=None, umat=None, vmat=None, transposedprj=False)

Apply a projection via solving a saddle point problem.

The theory is as follows. Consider the projection

\[P = I - M^{-1}J_1^T(J_1^TM^{-1}J_2)^{-1}J_2.\]

Then \(Pv\) can be obtained via

\[\begin{split}A^{-1}\begin{bmatrix} Pv \\ * \end{bmatrix} = \begin{bmatrix} Mv \\ 0 \end{bmatrix},\end{split}\]

where

\[\begin{split}A := \begin{bmatrix} M & J_1^T \\ J_2 & 0 \end{bmatrix}.\end{split}\]

And \(P^Tv\) can be obtained via

\[\begin{split}A^{-T}\begin{bmatrix} M^{-T}P^Tv \\ * \end{bmatrix} = \begin{bmatrix} v \\ 0 \end{bmatrix}.\end{split}\]

Parameters:

amat : (N,N) sparse matrix

left upper entry in the saddlepoint matrix

jmat : (M,N) sparse matrix

left lower entry

jmatT : (N,M) sparse matrix, optional

right upper entry, defaults to jmat.T

rhsv : (N,K) ndarray

array to be projected

umat, vmat : (N,L), (L,N) ndarrays or sparse matrices, optional

factored contribution to amat, default to None

transposedprj : boolean

whether to apply the transpose of the projection, defaults to False

Returns:

, : (N,K) ndarray

projected rhsv

lin_alg_utils.apply_sqrt_fromright(M, rhsa, output=None)

apply the sqrt of a mass matrix or other spd

TODO: cases for dense and sparse INPUTS

Parameters:

M : (N,N) sparse matrix

symmetric strictly positive definite

rhsa : (K,N) ndarray array or sparse matrix

array the inverse of M is to be applied to

output : string, optional

set to ‘sparse’ if rhsa has many zero rows to get the output as a sparse matrix

Returns:

, : (N,K) ndarray or sparse matrix

the sqrt of the inverse of M applied to rhsa from the left

lin_alg_utils.apply_invsqrt_fromright(M, rhsa, output=None)

apply the sqrt of the inverse of a mass matrix or other spd

TODO: cases for dense and sparse INPUTS

Parameters:

M : (N,N) sparse matrix

symmetric strictly positive definite

rhsa : (K,N) ndarray array or sparse matrix

array the inverse of M is to be applied to

output : string, optional

set to ‘sparse’ if rhsa has many zero rows to get the output as a sparse matrix

Returns:

, : (N,K) ndarray or sparse matrix

the sqrt of the inverse of M applied to rhsa from the left

lin_alg_utils.get_Sinv_smw(amat_lu, umat=None, vmat=None)

compute inverse of I-V*Ainv*U as it is needed for

the application of the Sherman-Morrison-Woodbury formula

Parameters:

amat_lu : callable f(v) or (N,N) sparse matrix

the main part of the matrix A-UV possibly lu-factored

umat, vmat : (N,L), (L,N) ndarrays or sparse matrices

factored contribution to amat_lu

Returns:

, : (L,L) ndarray

small inverse in the smw update

lin_alg_utils.solve_sadpnt_smw(amat=None, jmat=None, rhsv=None, jmatT=None, umat=None, vmat=None, rhsp=None, sadlu=None, return_alu=False, krylov=None, krpslvprms={}, krplsprms={})

solve a saddle point system

A - np.dot(U,V) J.T * X = rhsv J 0 rhsp by sparse direct solves

Parameters:

amat : (N,N) sparse matrix

left upper entry in the saddlepoint matrix

jmat : (M,N) sparse matrix

left lower entry

jmatT : (N,M) sparse matrix, optional

right upper entry, defaults to jmat.T

rhsv : (N,K) ndarray

upper part of right hand side

rhsp : (M,K) ndarray, optional

lower part of right hand side, defaults to zero

umat, vmat : (N,L), (L,N) ndarrays or sparse matrices, optional

factored contribution to amat, default to None

sadlu : callable f(v), optional

returns the inverse of the sadpoint matrix applied to v, defaults to None

return_alu : boolean, optional

whether to return the lu factored sadpoint matrix

Returns:

, : (N,K) ndarray

projected rhsv

, : f(v) callable, optional

lu decomposition of the saddlepoint matrix

lin_alg_utils.app_luinv_to_spmat(alu_solve, Z)

compute A.-1*Z where A comes factored

and with a solve routine for possibly complex Z

Parameters:

alu_solve : callable f(v)

returning a matrix inverse applied to v

Z : (N,K) ndarray, real or complex

the inverse is to be applied to

Returns:

, : (N,K) ndarray

matrix inverse applied to ndarray

lin_alg_utils.comp_sqfnrm_factrd_diff(zone, ztwo, ret_sing_norms=False)

compute the squared Frobenius norm of z1*z1.T - z2*z2.T

using the linearity traces and that tr(z1.dot(z2)) = tr(z2.dot(z1)) and that tr(z1.dot(z1.T)) is faster computed via (z1*z1.sum(-1)).sum()

lin_alg_utils.comp_sqfnrm_factrd_lyap_res(A, B, C)

compute the squared Frobenius norm of A*B.T + B*A.T + C*C.T

using the linearity traces and that tr.(z1.dot(z2)) = tr(z2.dot(z1)) and that tr(z1.dot(z1.T)) is faster computed via (z1*z1.sum(-1)).sum()

lin_alg_utils.comp_sqfnrm_factrd_sum(zone, ztwo, ret_sing_norms=False)

compute the squared Frobenius norm of z1*z1.T + z2*z2.T

using the linearity traces and that tr.(z1.dot(z2)) = tr(z2.dot(z1)) and that tr(z1.dot(z1.T)) is faster computed via (z1*z1.sum(-1)).sum()

lin_alg_utils.mm_dnssps(A, v)

compute A*v for sparse or dense A