Modelica.Math.Matrices.Utilities

Utility functions that should not be directly utilized by the user

Information

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).

Package Content

Name Description
Modelica.Math.Matrices.Utilities.continuousRiccatiIterative continuousRiccatiIterative Newton's method with exact line search for iterative solving continuous algebraic Riccati equation
Modelica.Math.Matrices.Utilities.discreteRiccatiIterative discreteRiccatiIterative Newton's method with exact line search for solving discrete algebraic Riccati equation
Modelica.Math.Matrices.Utilities.householderReflection householderReflection Reflect each of the vectors a_i of matrix A=[a_1, a_2, ..., a_n] on a plane with orthogonal vector u
Modelica.Math.Matrices.Utilities.householderSimilarityTransformation householderSimilarityTransformation Perform the similarity transformation S*A*S of matrix A with symmetric householder matrix S = I - 2u*u'
Modelica.Math.Matrices.Utilities.toUpperHessenberg toUpperHessenberg Transform a real square matrix A to upper Hessenberg form H by orthogonal similarity transformation: Q' * A * Q = H
Modelica.Math.Matrices.Utilities.eigenvaluesHessenberg eigenvaluesHessenberg Compute eigenvalues of an upper Hessenberg form matrix
Modelica.Math.Matrices.Utilities.reorderRSF reorderRSF Reorders a real Schur form to clusters of stable and unstable eigenvalues
Modelica.Math.Matrices.Utilities.findLocal_tk findLocal_tk Find a local minimizer tk to define the length of the step tk*Nk in continuousRiccatiIterative and discreteRiccatiIterative

Modelica.Math.Matrices.Utilities.continuousRiccatiIterative Modelica.Math.Matrices.Utilities.continuousRiccatiIterative

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

Information

Syntax

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

Description

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 X0 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].

References

  [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.

Example

     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;


See also

Matrices.Utilities.discreteRiccatiIterative
Matrices.continuousRiccati

Extends from Modelica.Icons.Function (Icon for functions).

Inputs

NameDescription
A[:, size(A, 1)]Matrix A of Riccati equation X*A + A'*X -X*G*X +Q = 0
B[size(A, 1), :]Matrix B in G = B*inv(R)*B'
R[size(B, 2), size(B, 2)]Matrix R in G = B*inv(R)*B'
Q[size(A, 1), size(A, 2)]Matrix Q of Riccati equation X*A + A'*X -X*G*X +Q = 0
X0[size(A, 1), size(A, 2)]Initial approximate solution for X*A + A'*X -X*G*X +Q = 0
maxStepsMaximal number of iteration steps
epsTolerance for stop criterion

Outputs

NameDescription
X[size(X0, 1), size(X0, 2)]Solution X of Riccati equation X*A + A'*X -X*G*X +Q = 0
rNorm of X*A + A'*X - X*G*X + Q, zero for exact solution

Modelica.Math.Matrices.Utilities.discreteRiccatiIterative Modelica.Math.Matrices.Utilities.discreteRiccatiIterative

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

Information

Syntax

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

Description

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 X0 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].

References

  [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.

Example

     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

See also

Matrices.Utilities.continuousRiccatiIterative
Matrices.discreteRiccati

Extends from Modelica.Icons.Function (Icon for functions).

Inputs

NameDescription
A[:, size(A, 1)]Matrix A of discrete Riccati equation
B[size(A, 1), :]Matrix B of discrete Riccati equation
R[size(B, 2), size(B, 2)]Matrix R of discrete Riccati equation
Q[size(A, 1), size(A, 2)]Matrix Q of discrete Riccati equation
X0[size(A, 1), size(A, 2)]Initial approximate solution discrete Riccati equation
maxStepsMaximal number of iteration steps
epsTolerance for stop criterion

Outputs

NameDescription
X[size(X0, 1), size(X0, 2)] 
r 

Modelica.Math.Matrices.Utilities.householderReflection Modelica.Math.Matrices.Utilities.householderReflection

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

Information

Syntax

Matrices.householderReflection(A,u);

Description

This function computes the Householder reflection (transformation)

Ar = Q*A
with
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.

Example

// 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]

See also

Matrices.Utilities.housholderSimilarityTransformation,
Vectors.Utilities.householderReflection,
Vectors.Utilities.householderVector

Extends from Modelica.Icons.Function (Icon for functions).

Inputs

NameDescription
A[:, :]Rectangular matrix
u[size(A, 1)]Householder vector

Outputs

NameDescription
RA[size(A, 1), size(A, 2)]Reflexion of A

Modelica.Math.Matrices.Utilities.householderSimilarityTransformation Modelica.Math.Matrices.Utilities.householderSimilarityTransformation

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

Information

Syntax

  As = Matrices.householderSimilarityTransformation(A,u);

Description

This function computes the Householder similarity transformation

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

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

Example

// 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]

See also

Matrices.Utilities.householderReflection,
Vectors.Utilities.householderReflection,
Vectors.Utilities.householderVector

Extends from Modelica.Icons.Function (Icon for functions).

Inputs

NameDescription
A[:, size(A, 1)]Square matrix A
u[size(A, 1)]Householder vector

Outputs

NameDescription
SAS[size(A, 1), size(A, 1)]Transformation of matrix A

Modelica.Math.Matrices.Utilities.toUpperHessenberg Modelica.Math.Matrices.Utilities.toUpperHessenberg

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

Information

Syntax

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

Description

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.

Example

 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]

See also

Matrices.hessenberg

Extends from Modelica.Icons.Function (Icon for functions).

Inputs

NameDescription
A[:, size(A, 1)]Square matrix A
iloLowest index where the original matrix had been Hessenbergform
ihiHighest index where the original matrix had been Hessenbergform

Outputs

NameDescription
H[size(A, 1), size(A, 2)]Upper Hessenberg form
V[size(A, 1), size(A, 2)]V=[v1,v2,..vn-1,0] with vi are vectors which define the elementary reflectors
tau[max(0, size(A, 1) - 1)]Scalar factors of the elementary reflectors
infoInformation of successful function call

Modelica.Math.Matrices.Utilities.eigenvaluesHessenberg Modelica.Math.Matrices.Utilities.eigenvaluesHessenberg

Compute eigenvalues of an upper Hessenberg form matrix

Information

Syntax

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

Description

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

Example

     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}

See also

Matrices.eigenValues, Matrices.hessenberg

Extends from Modelica.Icons.Function (Icon for functions).

Inputs

NameDescription
H[:, size(H, 1)]Hessenberg matrix H

Outputs

NameDescription
ev[size(H, 1), 2]Eigenvalues
info 

Modelica.Math.Matrices.Utilities.reorderRSF Modelica.Math.Matrices.Utilities.reorderRSF

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

Information

Syntax

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

Description

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.

Example

  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).

Inputs

NameDescription
T[:, :]Real Schur form
Q[:, size(T, 2)]Schur vector Matrix
alphaReal[size(T, 1)]Real part of eigenvalue=alphaReal+i*alphaImag
alphaImag[size(T, 1)]Imaginary part of eigenvalue=alphaReal+i*alphaImag
iscontinuousTrue if the according system is continuous. False for discrete systems

Outputs

NameDescription
To[size(T, 1), size(T, 2)]Reordered Schur form
Qo[size(T, 1), size(T, 2)]Reordered Schur vector matrix
wr[size(T, 2)]Reordered eigenvalues, real part
wi[size(T, 2)]Reordered eigenvalues, imaginary part

Modelica.Math.Matrices.Utilities.findLocal_tk Modelica.Math.Matrices.Utilities.findLocal_tk

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

Information

Syntax

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

Description

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

References

  [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.

See also

Matrices.Utilities.continuousRiccatiIterative
Matrices.Utilities.discreteRiccatiIterative

Extends from Modelica.Icons.Function (Icon for functions).

Inputs

NameDescription
Rk[:, size(Rk, 1)] 
Vk[size(Rk, 1), size(Rk, 2)] 

Outputs

NameDescription
tk 
Automatically generated Thu Dec 19 17:20:25 2019.