Utility functions that should not be directly utilized by the user

This package contains utility functions that are utilized by higher level matrix functions. These functions are usually not useful for an end-user.

Extends from `Modelica.Icons.UtilitiesPackage`

(Icon for utility packages).

Name | Description |
---|---|

`continuousRiccatiIterative` | Newton's method with exact line search for iterative solving continuous algebraic Riccati equation |

`discreteRiccatiIterative` | Newton's method with exact line search for solving discrete algebraic Riccati equation |

`eigenvaluesHessenberg` | Compute eigenvalues of an upper Hessenberg form matrix |

`findLocal_tk` | Find a local minimizer tk to define the length of the step tk*Nk in continuousRiccatiIterative and discreteRiccatiIterative |

`householderReflection` | Reflect each of the vectors a_i of matrix A=[a_1, a_2, ..., a_n] on a plane with orthogonal vector u |

`householderSimilarityTransformation` | Perform the similarity transformation S*A*S of matrix A with symmetric householder matrix S = I - 2u*u' |

`reorderRSF` | Reorders a real Schur form to clusters of stable and unstable eigenvalues |

`toUpperHessenberg` | Transform a real square matrix A to upper Hessenberg form H by orthogonal similarity transformation: Q' * A * Q = H |

Newton's method with exact line search for iterative solving continuous algebraic Riccati equation

X = Matrices.Utilities.continuousRiccatiIterative(A, B, R, Q, X0); (X, r) = Matrices.Utilities.continuousRiccatiIterative(A, B, R, Q, X0, maxSteps, eps);

This function provides a Newton-like method for solving continuous algebraic Riccati equations (care). It utilizes Exact Line Search to improve the sometimes erratic
convergence of Newton's method. Exact line search in this case means, that at each iteration `i`

a Newton step **delta**_i

X_i+1 =X_i +delta_i

is taken in the direction to minimize the Frobenius norm of the residual

r = ||X_i+1*A+A'*X_i+1 -X_i+1*G*X_i+1 +Q||.

with

-1 G =B*R*B'

The inputs "maxSteps" and "eps" specify the termination of the iteration. The iteration is terminated if either
maxSteps iteration steps have been performed or the relative change **delta**_i/**X**_i became smaller than eps.

With an appropriate initial value **X**0 a sufficiently accurate solution might be reach within a few iteration steps. Although a Lyapunov equation
of order `n`

(n is the order of the Riccati equation) is to be solved at each iteration step, the algorithm might be faster
than a direct method like Matrices.continuousRiccati, since direct methods have to solve the 2*n-order Hamiltonian
system equation.

The algorithm is taken from [1] and [2].

[1] Benner, P., Byers, R. An Exact Line Search Method for Solving Generalized Continuous-Time Algebraic Riccati Equations IEEE Transactions On Automatic Control, Vol. 43, No. 1, pp. 101-107, 1998. [2] Datta, B.N. Numerical Methods for Linear Control Systems Elsevier Academic Press, 2004.

A=[0.0, 1.0, 0.0, 0.0; 0.0, -1.890, 3.900e-01, -5.530; 0.0, -3.400e-02, -2.980, 2.430; 3.400e-02, -1.100e-03, -9.900e-01, -2.100e-01]; B=[ 0.0, 0.0; 3.600e-01, -1.60; -9.500e-01, -3.200e-02; 3.000e-02, 0.0]; R=[1, 0; 0, 1]; Q=[2.313, 2.727, 6.880e-01, 2.300e-02; 2.727, 4.271, 1.148, 3.230e-01; 6.880e-01, 1.148, 3.130e-01, 1.020e-01; 2.300e-02, 3.230e-01, 1.020e-01, 8.300e-02]; X0=identity(4); (X,r) = Matrices.Utilities.continuousRiccatiIterative(A, B, R, Q, X0); // X = [1.3239, 0.9015, 0.5466, -1.7672; 0.9015, 0.9607, 0.4334, -1.1989; 0.5466, 0.4334, 0.4605, -1.3633; -1.7672, -1.1989, -1.3633, 4.4612] // r = 2.48809423389491E-015 (,r) = Matrices.Utilities.continuousRiccatiIterative(A, B, R, Q, X0,4); // r = 0.0004;

Matrices.continuousRiccati

Extends from `Modelica.Icons.Function`

(Icon for functions).

Type | Name | Description |
---|---|---|

`Real` | `A[:,size(A, 1)]` | Matrix A of Riccati equation X*A + A'*X -X*G*X +Q = 0 |

`Real` | `B[size(A, 1),:]` | Matrix B in G = B*inv(R)*B' |

`Real` | `R[size(B, 2),size(B, 2)]` | Matrix R in G = B*inv(R)*B' |

`Real` | `Q[size(A, 1),size(A, 2)]` | Matrix Q of Riccati equation X*A + A'*X -X*G*X +Q = 0 |

`Real` | `X0[size(A, 1),size(A, 2)]` | Initial approximate solution for X*A + A'*X -X*G*X +Q = 0 |

`Integer` | `maxSteps` | Maximal number of iteration steps |

`Real` | `eps` | Tolerance for stop criterion |

Type | Name | Description |
---|---|---|

`Real` | `X[size(X0, 1),size(X0, 2)]` | Solution X of Riccati equation X*A + A'*X -X*G*X +Q = 0 |

`Real` | `r` | Norm of X*A + A'*X - X*G*X + Q, zero for exact solution |

Newton's method with exact line search for solving discrete algebraic Riccati equation

X = Matrices.Utilities.discreteRiccatiIterative(A, B, R, Q, X0); (X, r) = Matrices.Utilities.discreteRiccatiIterative(A, B, R, Q, X0, maxSteps, eps);

This function provides a Newton-like method for solving discrete-time algebraic Riccati equations. It uses Exact Line Search to improve the sometimes erratic
convergence of Newton's method. Exact line search in this case means, that at each iteration `i`

a Newton step **delta**_i

X_i+1 =X_i +delta_i

is taken in the direction to minimize the Frobenius norm of the residual

r = ||A'X_i+1*A-X_i+1 -A'X_i+1*G_i*X_i+1*A+Q||

with

-1 G_i =B*(R+B'*X_i*B) *B'

Output `r`

is the norm of the residual of the last iteration.

The inputs "maxSteps" and "eps" specify the termination of the iteration. The iteration is terminated if either
maxSteps iteration steps have been performed or the relative change **delta**_i/**X**_i became smaller than eps.

With an appropriate initial value **X**0 a sufficiently accurate solution might be reach with a few iteration steps. Although a Lyapunov equation of
order `n`

(n is the order of the Riccati equation) is to be solved at each iteration step, the algorithm might be faster
than a direct method like Matrices.discreteRiccati, since direct methods have to solve the 2*n-order Hamiltonian
system equation.
The algorithm is taken from [1] and [2].

[1] Benner, P., Byers, R. An Exact Line Search Method for Solving Generalized Continuous-Time Algebraic Riccati Equations IEEE Transactions On Automatic Control, Vol. 43, No. 1, pp. 101-107, 1998. [2] Datta, B.N. Numerical Methods for Linear Control Systems Elsevier Academic Press, 2004.

A = [0.9970, 0.0000, 0.0000, 0.0000; 1.0000, 0.0000, 0.0000, 0.0000; 0.0000, 1.0000, 0.0000, 0.0000; 0.0000, 0.0000, 1.0000, 0.0000]; B = [0.0150; 0.0000; 0.0000; 0.0000]; R = [0.2500]; Q = [0, 0, 0, 0; 0, 0, 0, 0; 0, 0, 0, 0; 0, 0, 0, 1]; X0=identity(4); (X,r) = Matrices.Utilities.discreteRiccatiIterative(A, B, R, Q, X0); // X = [30.625, 0.0, 0.0, 0.0; 0.0, 1.0, 0.0, 0.0; 0.0, 0.0, 1.0, 0.0; 0.0, 0.0, 0.0, 1.0]; // r = 3.10862446895044E-015

Matrices.discreteRiccati

Extends from `Modelica.Icons.Function`

(Icon for functions).

Type | Name | Description |
---|---|---|

`Real` | `A[:,size(A, 1)]` | Matrix A of discrete Riccati equation |

`Real` | `B[size(A, 1),:]` | Matrix B of discrete Riccati equation |

`Real` | `R[size(B, 2),size(B, 2)]` | Matrix R of discrete Riccati equation |

`Real` | `Q[size(A, 1),size(A, 2)]` | Matrix Q of discrete Riccati equation |

`Real` | `X0[size(A, 1),size(A, 2)]` | Initial approximate solution discrete Riccati equation |

`Integer` | `maxSteps` | Maximal number of iteration steps |

`Real` | `eps` | Tolerance for stop criterion |

Type | Name | Description |
---|---|---|

`Real` | `X[size(X0, 1),size(X0, 2)]` | |

`Real` | `r` |

Reflect each of the vectors a_i of matrix A=[a_1, a_2, ..., a_n] on a plane with orthogonal vector u

Matrices.householderReflection(A,u);

This function computes the Householder reflection (transformation)

withAr=Q*A

Q=I-2*u*u'/(u'*u)

where **u** is Householder vector, i.e., the normal vector of the reflection plane.

Householder reflection is widely used in numerical linear algebra, e.g., to perform QR decompositions.

// First step of QR decomposition import Modelica.Math.Vectors.Utilities; Real A[3,3] = [1,2,3; 3,4,5; 2,1,4]; Real Ar[3,3]; Real u[:]; u=Utilities.householderVector(A[:,1],{1,0,0}); // u= {0.763, 0.646, 0} Ar=householderReflection(A,u); // Ar = [-6.0828, -5.2608, -4.4388; // 0.0, -1.1508, -2.3016; // 0.0, 2.0, 0.0]

Matrices.Utilities.housholderSimilarityTransformation,

Vectors.Utilities.householderReflection,

Vectors.Utilities.householderVector

Extends from `Modelica.Icons.Function`

(Icon for functions).

Type | Name | Description |
---|---|---|

`Real` | `A[:,:]` | Rectangular matrix |

`Real` | `u[size(A, 1)]` | Householder vector |

Type | Name | Description |
---|---|---|

`Real` | `RA[size(A, 1),size(A, 2)]` | Reflexion of A |

Perform the similarity transformation S*A*S of matrix A with symmetric householder matrix S = I - 2u*u'

As = Matrices.householderSimilarityTransformation(A,u);

This function computes the Householder similarity transformation

withAs=S*A*S

S=I-2*u*u'/(u'*u).

This transformation is widely used for transforming non-symmetric matrices to a Hessenberg form.

// First step of Hessenberg decomposition import Modelica.Math.Vectors.Utilities; Real A[4,4] = [1,2,3,4; 3,4,5,6; 9,8,7,6; 1,2,0,0]; Real Ar[4,4]; Real u[4]={0,0,0,0}; u[2:4]=Utilities.householderVector(A[2:4,1],{1,0,0}); // u= = {0, 0.8107, 0.5819, 0.0647} Ar=householderSimilarityTransformation(A,u); // Ar = [1.0, -3.8787, -1.2193, 3.531; -9.5394, 11.3407, 6.4336, -5.9243; 0.0, 3.1307, 0.7525, -3.3670; 0.0, 0.8021, -1.1656, -1.0932]

Matrices.Utilities.householderReflection,

Vectors.Utilities.householderReflection,

Vectors.Utilities.householderVector

Extends from `Modelica.Icons.Function`

(Icon for functions).

Type | Name | Description |
---|---|---|

`Real` | `A[:,size(A, 1)]` | Square matrix A |

`Real` | `u[size(A, 1)]` | Householder vector |

Type | Name | Description |
---|---|---|

`Real` | `SAS[size(A, 1),size(A, 1)]` | Transformation of matrix A |

Transform a real square matrix A to upper Hessenberg form H by orthogonal similarity transformation: Q' * A * Q = H

H = Matrices.Utilities.toUpperHessenberg(A); (H, V, tau, info) = Matrices.Utilities.toUpperHessenberg(A,ilo, ihi);

Function **toUpperHessenberg** computes a upper Hessenberg form **H** of a matrix **A** by orthogonal similarity transformation: **Q**' * **A** * **Q** = **H**.
With the optional inputs ilo and ihi, also partial transformation is possible. The function calls LAPACK function DGEHRD.
See Matrices.LAPACK.dgehrd for more information about the additional outputs V, tau, info and inputs ilo, ihi.

A = [1, 2, 3; 6, 5, 4; 1, 0, 0]; H = toUpperHessenberg(A); results in: H = [1.0, -2.466, 2.630; -6.083, 5.514, -3.081; 0.0, 0.919, -0.514]

Extends from `Modelica.Icons.Function`

(Icon for functions).

Type | Name | Description |
---|---|---|

`Real` | `A[:,size(A, 1)]` | Square matrix A |

`Integer` | `ilo` | Lowest index where the original matrix had been Hessenbergform |

`Integer` | `ihi` | Highest index where the original matrix had been Hessenbergform |

Type | Name | Description |
---|---|---|

`Real` | `H[size(A, 1),size(A, 2)]` | Upper Hessenberg form |

`Real` | `V[size(A, 1),size(A, 2)]` | V=[v1,v2,..vn-1,0] with vi are vectors which define the elementary reflectors |

`Real` | `tau[max(0, size(A, 1) - 1)]` | Scalar factors of the elementary reflectors |

`Integer` | `info` | Information of successful function call |

Compute eigenvalues of an upper Hessenberg form matrix

ev = Matrices.Utilities.eigenvaluesHessenberg(H); (ev, info) = Matrices.Utilities.eigenvaluesHessenberg(H);

This function computes the eigenvalues of a Hessenberg form matrix. Transformation to Hessenberg form is the first step in eigenvalue computation for arbitrary matrices with QR decomposition. This step can be skipped if the matrix has already Hessenberg form.

The function uses the LAPACK-routine dhseqr. Output `info`

is 0 for a successful call of this
function.

See Matrices.LAPACK.dhseqr for details

Real A[3,3] = [1,2,3; 9,8,7; 0,1,0]; Real ev[3,2]; ev := Matrices.Utilities.eigenvaluesHessenberg(A); // ev = [10.7538, 0.0; -0.8769, 1.0444; -0.8769, -1.0444] // = {10.7538, -0.8769 +- i*1.0444}

Matrices.eigenValues, Matrices.hessenberg

Extends from `Modelica.Icons.Function`

(Icon for functions).

Type | Name | Description |
---|---|---|

`Real` | `H[:,size(H, 1)]` | Hessenberg matrix H |

Type | Name | Description |
---|---|---|

`Real` | `ev[size(H, 1),2]` | Eigenvalues |

`Integer` | `info` |

Reorders a real Schur form to clusters of stable and unstable eigenvalues

To = Matrices.Utilities.reorderRSF(T, Q, alphaReal, alphaImag); (To, Qo, wr, wi) = Matrices.Utilities.reorderRSF(T, Q, alphaReal, alphaImag, iscontinuous);

Function **reorderRSF**() reorders a real Schur form such that the stable eigenvalues of
the system are in the 1-by-1 and 2-by-2 diagonal blocks of the block **upper** triangular matrix.
If the Schur form is referenced to a continuous system the staple eigenvalues are in the left complex half plane.
The stable eigenvalues of a discrete system are inside the complex unit circle.

This function is used for example to solve algebraic Riccati equations
(continuousRiccati,
discreteRiccati). In this context the Schur form
as well as the corresponding eigenvalues and the transformation matrix **Q** are known, why the eigenvalues and the transformation matrix are inputs to **reorderRSF()**.

The Schur vector matrix **Qo** is also reordered according to **To**. The vectors **wr** and **wi** contains the real and imaginary parts of the
reordered eigenvalues respectively.

T := [-1,2, 3,4; 0,2, 6,5; 0,0,-3,5; 0,0, 0,6]; To := Matrices.Utilities.reorderRSF(T,identity(4),{-1, 2, -3, 6},{0, 0, 0, 0}, true); // To = [-1.0, -0.384, 3.585, 4.0; // 0.0, -3.0, 6.0, 0.64; // 0.0, 0.0, 2.0, 7.04; // 0.0, 0.0, 0.0, 6.0]

See also Matrices.realSchur

Extends from `Modelica.Icons.Function`

(Icon for functions).

Type | Name | Description |
---|---|---|

`Real` | `T[:,:]` | Real Schur form |

`Real` | `Q[:,size(T, 2)]` | Schur vector Matrix |

`Real` | `alphaReal[size(T, 1)]` | Real part of eigenvalue=alphaReal+i*alphaImag |

`Real` | `alphaImag[size(T, 1)]` | Imaginary part of eigenvalue=alphaReal+i*alphaImag |

`Boolean` | `iscontinuous` | True if the according system is continuous. False for discrete systems |

Type | Name | Description |
---|---|---|

`Real` | `To[size(T, 1),size(T, 2)]` | Reordered Schur form |

`Real` | `Qo[size(T, 1),size(T, 2)]` | Reordered Schur vector matrix |

`Real` | `wr[size(T, 2)]` | Reordered eigenvalues, real part |

`Real` | `wi[size(T, 2)]` | Reordered eigenvalues, imaginary part |

Find a local minimizer tk to define the length of the step tk*Nk in continuousRiccatiIterative and discreteRiccatiIterative

tk = Matrices.Utilities.findLocal_tk(Rk, Vk);

Function `findLocal_tk()`

is an auxiliary function called in iterative solver for algebraic Riccati equation based on Newton's method with
exact line search like continuousRiccatiIterative

and discreteRiccatiIterative.

The function computes the local minimum of the function f_k(t_k)

f_k(t_k) = alpha_k*(1-t_k)^2 + 2*beta_k*(1-t)*t^2 + gamma_k*t^4

by calculating the zeros of the derivation d f_k/d t_k. It is known that the function f_k(t_k) has a local minimum at some value t_k_min in [0, 2].

With t_k_min the norm of the next residual of the algorithm will be minimized.

See [1] for more information

[1] Benner, P., Byers, R. An Exact Line Search Method for Solving Generalized Continuous-Time Algebraic Riccati Equations IEEE Transactions On Automatic Control, Vol. 43, No. 1, pp. 101-107, 1998.

Matrices.Utilities.discreteRiccatiIterative

Extends from `Modelica.Icons.Function`

(Icon for functions).

Type | Name | Description |
---|---|---|

`Real` | `Rk[:,size(Rk, 1)]` | |

`Real` | `Vk[size(Rk, 1),size(Rk, 2)]` |

Type | Name | Description |
---|---|---|

`Real` | `tk` |