Modelica_LinearSystems2.Math.Matrices

Additional functions for Modelica.Math

Information

Extends from Modelica.Icons.Library (Icon for library).

Package Content

NameDescription
Modelica_LinearSystems2.Math.Matrices.LAPACK LAPACK  
Modelica_LinearSystems2.Math.Matrices.Examples Examples  
Modelica_LinearSystems2.Math.Matrices.care care Solution of continuous-time algebraic Riccati equations
Modelica_LinearSystems2.Math.Matrices.dare dare Solution of discrete-time algebraic Riccati equations
Modelica_LinearSystems2.Math.Matrices.det det Determinant of a matrix (computed by LU decomposition)
Modelica_LinearSystems2.Math.Matrices.eigenValues eigenValues Compute eigenvalues and eigenvectors for a real, nonsymmetric matrix
Modelica_LinearSystems2.Math.Matrices.equalityLeastSquares equalityLeastSquares Solve a linear equality constrained least squares problem
Modelica_LinearSystems2.Math.Matrices.fliplr fliplr flip the columns of a matrix in left/right direction
Modelica_LinearSystems2.Math.Matrices.flipud flipud flip the columns of a matrix in up/down direction
Modelica_LinearSystems2.Math.Matrices.frobeniusNorm frobeniusNorm Return the Frobenius norm of a matrix
Modelica_LinearSystems2.Math.Matrices.fromFile fromFile Read matrix from a matlab file
Modelica_LinearSystems2.Math.Matrices.generalizedEigenvaluesTriangular generalizedEigenvaluesTriangular Compute invariant zeros of linear state space system with a genralized system matrix [A, B, C, D] which is of upper Hessenberg form
Modelica_LinearSystems2.Math.Matrices.hessenberg hessenberg Compute an upper Hessenberg matrix by repeatedly applicated householder similarity transformation
Modelica_LinearSystems2.Math.Matrices.householderReflexion householderReflexion reflect each of the vectors ai of matrix A=[a1, a2, ..., an] on a plane with orthogonal vector u
Modelica_LinearSystems2.Math.Matrices.householderSimilarityTransformation householderSimilarityTransformation Calculate the similarity transformation SAS of matrix A with householder matrix S = I - 2u*u'
Modelica_LinearSystems2.Math.Matrices.leastSquares leastSquares Solve overdetermined or underdetermined real system of linear equations A*x=b in a least squares sense (A may be rank deficient)
Modelica_LinearSystems2.Math.Matrices.LU LU LU decomposition of square or rectangular matrix
Modelica_LinearSystems2.Math.Matrices.LU_solve LU_solve Solve real system of linear equations P*L*U*x=b with a b vector and an LU decomposition (from LU(..))
Modelica_LinearSystems2.Math.Matrices.LU_solve2 LU_solve2 Solve real system of linear equations P*L*U*X=B with a B vector and an LU decomposition (from LU(..))
Modelica_LinearSystems2.Math.Matrices.lyapunov lyapunov Solution of continuous-time Lyapunov equation X*A + A'*X = C
Modelica_LinearSystems2.Math.Matrices.norm norm Returns the norm of a matrix
Modelica_LinearSystems2.Math.Matrices.orthogonalQ orthogonalQ generates a real orthogonal matrix Q which is defined as the product of IHI-ILO elementary reflectors of order N, as returned by DGEHRD
Modelica_LinearSystems2.Math.Matrices.printMatrix printMatrix print matrix
Modelica_LinearSystems2.Math.Matrices.QR QR QR decomposition of a square matrix without column pivoting (A = Q*R)
Modelica_LinearSystems2.Math.Matrices.rsf rsf Computes the real Schur form (RSF) of a square matrix
Modelica_LinearSystems2.Math.Matrices.solve solve Solve real system of linear equations A*x=b with a b vector (Gaussian elemination with partial pivoting)
Modelica_LinearSystems2.Math.Matrices.solve2 solve2 Solve real system of linear equations A*X=B with a B matrix (Gaussian elemination with partial pivoting)
Modelica_LinearSystems2.Math.Matrices.toUpperHessenberg toUpperHessenberg transform a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation: Q' * A * Q = H
Modelica_LinearSystems2.Math.Matrices.trace trace tarce(A) is the sum of the diagonal elements of A

Modelica_LinearSystems2.Math.Matrices.care

Solution of continuous-time algebraic Riccati equations

Information


 
 
Function care computes the solution X of the continuous-time algebraic Riccati equation
 Q + A'*X + X*A - X*G*X = 0 
with
       -1
G = B*R *B' 
using the Schur vector approach proposed by Laub [1].

It is assumed that Q is symmetric and positve semidefinite and R is symmetric, nonsingular and positive definite, (A,B) is stabilizable and (A,Q) is detectable.

The assumptions are not checked in this function

The assumptions guarantee that Hamiltonian matrix

H = [A, -G; -Q, -A']
has no pure imaginary eigenvalue and can be put to an ordered real Schur form
U'*H*U = S = [S11, S12; 0, S22]
with orthogonal similarity transformation U. S is ordered in such a way, that S11 contains the n stable eigenvalues of the closed loop system with system matrix
       -1
A - B*R *B'*X
If U is partitioned to
U = [U11, U12; U21, U22]
with dimenstions according to S, the solution X can be calculated by
X*U11 = U21.
The algorithm uses LAPACK routines dgehrd (to compute the upper Hessenberg matrix of H), dorghr (to calculate the orthogonal matrix from the elementary reflectors as returned from dgehrd), dhseqr (to put transformed H to Schur form and to calculate the eigenvalues of the closed loop system) and dtrsen (to compute the ordered real Schur form and matrix U).

References

  [1] Laub, A.J.
      A Schur Method for Solving Algebraic Riccati equations.
      IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

Inputs

NameDescription
A[:, size(A, 1)] 
B[size(A, 1), :] 
R[size(B, 2), size(B, 2)] 
Q[size(A, 1), size(A, 1)] 
refine 

Outputs

NameDescription
X[size(A, 1), size(A, 2)]stabilizing solution of CARE
ev[size(A, 1)]eigenvalues of the closed loop system

Modelica_LinearSystems2.Math.Matrices.dare

Solution of discrete-time algebraic Riccati equations

Information


Function dare computes the solution X of the discrete-time algebraic Riccati equation
                                 -1
 X = A'*X*A - A'*X*B*(R + B'*X*B)  *B'*X*A + Q
using the Schur vector approach proposed by Laub [1].

It is assumed that Q is symmetric and positve semidefinite and R is symmetric, nonsingular and positive definite, (A,B) is stabilizable and (A,Q) is detectable.

The assumptions are not checked in this function

The assumptions guarantee that the Hamiltonian matrix

      -1  -1       -1        -1
H = [A, -A  *G; Q*A, A' + Q*A  *G ]
with
       -1
G = B*R *B' 
has no eigenvalue on the unit circle and can be put to an ordered real Schur form
U'*H*U = S = [S11, S12; 0, S22]
with orthogonal similarity transformation U. S is ordered in such a way, that S11 contains the n stable eigenvalues of the closed loop system with system matrix
                  -1
A - B*(R + B'*X*B)  *B'*X*A
If U is partitioned to
U = [U11, U12; U21, U22]
according to S, the solution X can be calculated by
X*U11 = U21.
The algorithm uses LAPACK routines dgehrd (to compute the upper Hessenberg matrix of H), dorghr (to calculate the orthogonal matrix from the elementary reflectors as returned from dgehrd), dhseqr (to put transformed H to Schur form and to calculate the eigenvalues of the closed loop system) and dtrsen (to compute the ordered real Schur form and matrix U).

References

  [1] Laub, A.J.
      A Schur Method for Solving Algebraic Riccati equations.
      IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

Inputs

NameDescription
A[:, size(A, 1)] 
B[size(A, 1), :] 
R[size(B, 2), size(B, 2)] 
Q[size(A, 1), size(A, 1)] 

Outputs

NameDescription
S[size(A, 1), size(A, 2)]orthogonal matrix of the Schur vectors associated to ordered rsf
ev[size(A, 1)]eigenvalues of the closed loop system

Modelica_LinearSystems2.Math.Matrices.det Modelica_LinearSystems2.Math.Matrices.det

Determinant of a matrix (computed by LU decomposition)

Information


Syntax

Matrices.det(A);

Description

This function call returns the determinant of matrix A computed by a LU decomposition. Usally, this function should never be used, because there are nearly always better numerical algorithms as by computing the determinant. E.g., use function Matrices.rank to compute the rank of a matrix.

See also

Matrices.rank, Matrices.solve

Extends from Modelica.Icons.Function (Icon for a function).

Inputs

NameDescription
A[:, size(A, 1)] 

Outputs

NameDescription
resultDeterminant of matrix A

Modelica_LinearSystems2.Math.Matrices.eigenValues Modelica_LinearSystems2.Math.Matrices.eigenValues

Compute eigenvalues and eigenvectors for a real, nonsymmetric matrix

Information


Syntax

                eigenvalues = Matrices.eigenValues(A);
(eigenvalues, eigenvectors) = Matrices.eigenValues(A);

Description

This function call returns the eigenvalues and optionally the (right) eigenvectors of a square matrix A. The first column of "eigenvalues" contains the real and the second column contains the imaginary part of the eigenvalues. If the i-th eigenvalue has no imaginary part, then eigenvectors[:,i] is the corresponding real eigenvector. If the i-th eigenvalue has an imaginary part, then eigenvalues[i+1,:] is the conjugate complex eigenvalue and eigenvectors[:,i] is the real and eigenvectors[:,i+1] is the imaginary part of the eigenvector of the i-th eigenvalue. With function Matrices.eigenValueMatrix, a real block diagonal matrix is constructed from the eigenvalues such that

A = eigenvectors * eigenValueMatrix(eigenvalues) * inv(eigenvectors)

provided the eigenvector matrix "eigenvectors" can be inverted (an inversion is possible, if all eigenvalues are different and no eigenvalue is zero).

Example

  Real A[3,3] = [1,2,3; 
                 3,4,5;
                 2,1,4];
  Real eval;
algorithm
  eval := Matrices.eigenValues(A);  // eval = [-0.618, 0; 
                                    //          8.0  , 0;
                                    //          1.618, 0];

i.e., matrix A has the 3 real eigenvalues -0.618, 8, 1.618.

See also

Matrices.eigenValueMatrix, Matrices.singularValues

Extends from Modelica.Icons.Function (Icon for a function).

Inputs

NameDescription
A[:, size(A, 1)]Matrix

Outputs

NameDescription
eigenvalues[size(A, 1), 2]Eigenvalues of matrix A (Re: first column, Im: second column)
leftEigenvectors[size(A, 1), size(A, 2)]Real-valued eigenvector matrix
rightEigenvectors[size(A, 1), size(A, 2)]Real-valued eigenvector matrix

Modelica_LinearSystems2.Math.Matrices.equalityLeastSquares Modelica_LinearSystems2.Math.Matrices.equalityLeastSquares

Solve a linear equality constrained least squares problem

Information


Syntax

x = Matrices.equalityLeastSquares(A,a,B,b);

Description

This function returns the solution x of the linear equality-constrained least squares problem:

min|A*x - a|^2 over x, subject to B*x = b

It is required that the dimensions of A and B fulfill the following relationship:

size(B,1) ≤ size(A,2) ≤ size(A,1) + size(B,1)

Note, the solution is computed with the LAPACK function "dgglse" using the generalized RQ factorization under the assumptions that B has full row rank (= size(B,1)) and the matrix [A;B] has full column rank (= size(A,2)). In this case, the problem has a unique solution.

Extends from Modelica.Icons.Function (Icon for a function).

Inputs

NameDescription
A[:, :]Minimize |A*x - a|^2
a[size(A, 1)] 
B[:, size(A, 2)]subject to B*x=b
b[size(B, 1)] 

Outputs

NameDescription
x[size(A, 2)]solution vector

Modelica_LinearSystems2.Math.Matrices.fliplr

flip the columns of a matrix in left/right direction

Inputs

NameDescription
A[:, :]Matrix to be fliped

Outputs

NameDescription
Aflip[size(A, 1), size(A, 2)]fliped matrix

Modelica_LinearSystems2.Math.Matrices.flipud

flip the columns of a matrix in up/down direction

Inputs

NameDescription
A[:, :]Matrix to be fliped

Outputs

NameDescription
Aflip[size(A, 1), size(A, 2)]fliped matrix

Modelica_LinearSystems2.Math.Matrices.frobeniusNorm Modelica_LinearSystems2.Math.Matrices.frobeniusNorm

Return the Frobenius norm of a matrix

Information

Extends from Modelica.Icons.Function (Icon for a function).

Inputs

NameDescription
A[:, :]Input matrix

Outputs

NameDescription
resultfrobenius norm of matrix A

Modelica_LinearSystems2.Math.Matrices.fromFile

Read matrix from a matlab file

Inputs

NameDescription
fileName 
matrixNameName of the matrix

Outputs

NameDescription
A[n, m] 

Modelica_LinearSystems2.Math.Matrices.generalizedEigenvaluesTriangular

Compute invariant zeros of linear state space system with a genralized system matrix [A, B, C, D] which is of upper Hessenberg form

Information

This function is an interface to LAPACK routine DHGEQZ to calculate invariant
zeros of systems with generalized system matrices of upper Hessenberg form.
DHGEQZ is described below:
 
 
 
     Purpose  
   ==========================================================
 
   DHGEQZ implements a single-/double-shift version of the QZ method for  
   finding the generalized eigenvalues  
 
   w(j)=(ALPHAR(j) + i*ALPHAI(j))/BETAR(j)   of the equation  
 
        det( A - w(i) B ) = 0  
 
   In addition, the pair A,B may be reduced to generalized Schur form:  
   B is upper triangular, and A is block upper triangular, where the  
   diagonal blocks are either 1-by-1 or 2-by-2, the 2-by-2 blocks having  
   complex generalized eigenvalues (see the description of the argument  
   JOB.)  
 
   If JOB='S', then the pair (A,B) is simultaneously reduced to Schur  
   form by applying one orthogonal tranformation (usually called Q) on  
   the left and another (usually called Z) on the right.  The 2-by-2  
   upper-triangular diagonal blocks of B corresponding to 2-by-2 blocks  
   of A will be reduced to positive diagonal matrices.  (I.e.,  
   if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and  
   B(j+1,j+1) will be positive.)  
 
   If JOB='E', then at each iteration, the same transformations  
   are computed, but they are only applied to those parts of A and B  
   which are needed to compute ALPHAR, ALPHAI, and BETAR.  
 
   If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal  
   transformations used to reduce (A,B) are accumulated into the arrays  
   Q and Z s.t.:  
 
        Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)*  
        Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)*  
 
   Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix  
        Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),  
        pp. 241--256.  
 
   Arguments  
   =========  
 
   JOB     (input) CHARACTER*1  
           = 'E': compute only ALPHAR, ALPHAI, and BETA.  A and B will  
                  not necessarily be put into generalized Schur form.  
           = 'S': put A and B into generalized Schur form, as well  
                  as computing ALPHAR, ALPHAI, and BETA.  
 
   COMPQ   (input) CHARACTER*1  
           = 'N': do not modify Q.  
           = 'V': multiply the array Q on the right by the transpose of  
                  the orthogonal tranformation that is applied to the  
                  left side of A and B to reduce them to Schur form.  
           = 'I': like COMPQ='V', except that Q will be initialized to  
                  the identity first.  
 
   COMPZ   (input) CHARACTER*1  
           = 'N': do not modify Z.  
           = 'V': multiply the array Z on the right by the orthogonal  
                  tranformation that is applied to the right side of  
                  A and B to reduce them to Schur form.  
           = 'I': like COMPZ='V', except that Z will be initialized to  
                  the identity first.  
 
   N       (input) INTEGER  
           The order of the matrices A, B, Q, and Z.  N >= 0.  
 
   ILO     (input) INTEGER  
   IHI     (input) INTEGER  
           It is assumed that A is already upper triangular in rows and  
           columns 1:ILO-1 and IHI+1:N.  
           1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.  
 
   A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)  
           On entry, the N-by-N upper Hessenberg matrix A.  Elements  
           below the subdiagonal must be zero.  
           If JOB='S', then on exit A and B will have been  
              simultaneously reduced to generalized Schur form.  
           If JOB='E', then on exit A will have been destroyed.  
              The diagonal blocks will be correct, but the off-diagonal  
              portion will be meaningless.  
 
   LDA     (input) INTEGER  
           The leading dimension of the array A.  LDA >= max( 1, N ).  
 
   B       (input/output) DOUBLE PRECISION array, dimension (LDB, N)  
           On entry, the N-by-N upper triangular matrix B.  Elements  
           below the diagonal must be zero.  2-by-2 blocks in B  
           corresponding to 2-by-2 blocks in A will be reduced to  
           positive diagonal form.  (I.e., if A(j+1,j) is non-zero,  
           then B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be  
           positive.)  
           If JOB='S', then on exit A and B will have been  
              simultaneously reduced to Schur form.  
           If JOB='E', then on exit B will have been destroyed.  
              Elements corresponding to diagonal blocks of A will be  
              correct, but the off-diagonal portion will be meaningless.  
 
   LDB     (input) INTEGER  
           The leading dimension of the array B.  LDB >= max( 1, N ).  
 
   ALPHAR  (output) DOUBLE PRECISION array, dimension (N)  
           ALPHAR(1:N) will be set to real parts of the diagonal  
           elements of A that would result from reducing A and B to  
           Schur form and then further reducing them both to triangular  
           form using unitary transformations s.t. the diagonal of B  
           was non-negative real.  Thus, if A(j,j) is in a 1-by-1 block  
           (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=A(j,j).  
           Note that the (real or complex) values  
           (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the  
           generalized eigenvalues of the matrix pencil A - wB.  
 
   ALPHAI  (output) DOUBLE PRECISION array, dimension (N)  
           ALPHAI(1:N) will be set to imaginary parts of the diagonal  
           elements of A that would result from reducing A and B to  
           Schur form and then further reducing them both to triangular  
           form using unitary transformations s.t. the diagonal of B  
           was non-negative real.  Thus, if A(j,j) is in a 1-by-1 block  
           (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=0.  
           Note that the (real or complex) values  
           (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the  
           generalized eigenvalues of the matrix pencil A - wB.  
 
   BETA    (output) DOUBLE PRECISION array, dimension (N)  
           BETA(1:N) will be set to the (real) diagonal elements of B  
           that would result from reducing A and B to Schur form and  
           then further reducing them both to triangular form using  
           unitary transformations s.t. the diagonal of B was  
           non-negative real.  Thus, if A(j,j) is in a 1-by-1 block  
           (i.e., A(j+1,j)=A(j,j+1)=0), then BETA(j)=B(j,j).  
           Note that the (real or complex) values  
           (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the  
           generalized eigenvalues of the matrix pencil A - wB.  
           (Note that BETA(1:N) will always be non-negative, and no  
           BETAI is necessary.)  
 
   Q       (input/output) DOUBLE PRECISION array, dimension (LDQ, N)  
           If COMPQ='N', then Q will not be referenced.  
           If COMPQ='V' or 'I', then the transpose of the orthogonal  
              transformations which are applied to A and B on the left  
              will be applied to the array Q on the right.  
 
   LDQ     (input) INTEGER  
           The leading dimension of the array Q.  LDQ >= 1.  
           If COMPQ='V' or 'I', then LDQ >= N.  
 
   Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)  
           If COMPZ='N', then Z will not be referenced.  
           If COMPZ='V' or 'I', then the orthogonal transformations  
              which are applied to A and B on the right will be applied  
              to the array Z on the right.  
 
   LDZ     (input) INTEGER  
           The leading dimension of the array Z.  LDZ >= 1.  
           If COMPZ='V' or 'I', then LDZ >= N.  
 
   WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)  
           On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.  
 
   LWORK   (input) INTEGER  
           The dimension of the array WORK.  LWORK >= max(1,N).  
 
           If LWORK = -1, then a workspace query is assumed; the routine  
           only calculates the optimal size of the WORK array, returns  
           this value as the first entry of the WORK array, and no error  
           message related to LWORK is issued by XERBLA.  
 
   INFO    (output) INTEGER  
           = 0: successful exit  
           < 0: if INFO = -i, the i-th argument had an illegal value  
           = 1,...,N: the QZ iteration did not converge.  (A,B) is not  
                      in Schur form, but ALPHAR(i), ALPHAI(i), and  
                      BETA(i), i=INFO+1,...,N should be correct.  
           = N+1,...,2*N: the shift calculation failed.  (A,B) is not  
                      in Schur form, but ALPHAR(i), ALPHAI(i), and  
                      BETA(i), i=INFO-N+1,...,N should be correct.  
           > 2*N:     various "impossible" errors.  
 
   Further Details  
   ===============  
 
   Iteration counters:  
 
   JITER  -- counts iterations.  
   IITER  -- counts iterations run since ILAST was last  
             changed.  This is therefore reset only when a 1-by-1 or  
             2-by-2 block deflates off the bottom.  
 
   =====================================================================  

Inputs

NameDescription
A[:, size(A, 1)] 
B[size(A, 1), size(A, 1)] 

Outputs

NameDescription
alphaReal[size(A, 1)]Real part of alpha (eigenvalue=(alphaReal+i*alphaImag)/beta)
alphaImag[size(A, 1)]Imaginary part of alpha
beta[size(A, 1)]Denominator of eigenvalue
info 

Modelica_LinearSystems2.Math.Matrices.hessenberg

Compute an upper Hessenberg matrix by repeatedly applicated householder similarity transformation

Information


 
This function computes the Hessenberg matrix of matrix A by repetitive application of Householder similarity transformation 
 
    Ai+1 = (I-2*u_i*u_i')*Ai*(I-2*u_i*u_i')
with Householder vector u_i.

The elementary transformations can be subsumed under

 A -> Q*A*Q
and Q*A*Q is Hessenberg matrix.

Modelica_LinearSystems2.Math.Matrices.hess uses LAPACK routine dgehrd. In contrast to this function Modelica_LinearSystems2.Math.Matrices.Internal.hessenberg does not use any LAPACK routine.

Inputs

NameDescription
H[size(H, 1), size(H, 2)] 

Outputs

NameDescription
Ht[size(H, 1), size(H, 2)] 

Modelica_LinearSystems2.Math.Matrices.householderReflexion

reflect each of the vectors ai of matrix A=[a1, a2, ..., an] on a plane with orthogonal vector u

Inputs

NameDescription
A[:, :] 
u[size(A, 1)]householder vector

Outputs

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

Modelica_LinearSystems2.Math.Matrices.householderSimilarityTransformation

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

Inputs

NameDescription
A[:, size(A, 1)] 
u[size(A, 1)]householder vector

Outputs

NameDescription
SAS[size(A, 1), size(A, 1)] 

Modelica_LinearSystems2.Math.Matrices.leastSquares Modelica_LinearSystems2.Math.Matrices.leastSquares

Solve overdetermined or underdetermined real system of linear equations A*x=b in a least squares sense (A may be rank deficient)

Information


Syntax

x = Matrices.leastSquares(A,b);

Description

A linear system of equations A*x = b has no solutions or infinitely many solutions if A is not square. Function "leastSquares" returns a solution in a least squarse sense:

  size(A,1) > size(A,2):  returns x such that |A*x - b|^2 is a minimum
  size(A,1) = size(A,2):  returns x such that A*x = b
  size(A,1) < size(A,2):  returns x such that |x|^2 is a minimum for all 
                          vectors x that fulfill A*x = b

Note, the solution is computed with the LAPACK function "dgelsx", i.e., QR or LQ factorization of A with column pivoting. If A does not have full rank, the solution is not unique and from the infinitely many solutions the one is selected that minimizes both |x|^2 and |A*x - b|^2.

Extends from Modelica.Icons.Function (Icon for a function).

Inputs

NameDescription
A[:, :]Matrix A
b[size(A, 1)]Vector b

Outputs

NameDescription
x[size(A, 2)]Vector x such that min|A*x-b|^2 if size(A,1) >= size(A,2) or min|x|^2 and A*x=b, if size(A,1) < size(A,2)

Modelica_LinearSystems2.Math.Matrices.LU Modelica_LinearSystems2.Math.Matrices.LU

LU decomposition of square or rectangular matrix

Information


Syntax

(LU, pivots)       = Matrices.LU(A);
(LU, pivots, info) = Matrices.LU(A);

Description

This function call returns the LU decomposition of a "Real[m,n]" matrix A, i.e.,

P*L*U = A

where P is a permutation matrix (implicitely defined by vector pivots), L is a lower triangular matrix with unit diagonal elements (lower trapezoidal if m > n), and U is an upper triangular matrix (upper trapezoidal if m < n). Matrices L and U are stored in the returned matrix LU (the diagonal of L is not stored). With the companion function Matrices.LU_solve, this decomposition can be used to solve linear systems (P*L*U)*x = b with different right hand side vectors b. If a linear system of equations with just one right hand side vector b shall be solved, it is more convenient to just use the function Matrices.solve.

The optional third (Integer) output argument has the following meaning: successful exit
info = 0:
info > 0: if info = i, U[i,i] is exactly zero. The factorization has been completed,
but the factor U is exactly singular, and division by zero will occur
if it is used to solve a system of equations.

The LU factorization is computed with the LAPACK function "dgetrf", i.e., by Gaussian elemination using partial pivoting with row interchanges. Vector "pivots" are the pivot indices, i.e., for 1 ≤ i ≤ min(m,n), row i of matrix A was interchanged with row pivots[i].

Example

  Real A[3,3] = [1,2,3; 
                 3,4,5;
                 2,1,4];
  Real b1[3] = {10,22,12};
  Real b2[3] = { 7,13,10};
  Real    LU[3,3];
  Integer pivots[3];
  Real    x1[3];
  Real    x2[3];
algorithm
  (LU, pivots) := Matrices.LU(A);
  x1 := Matrices.LU_solve(LU, pivots, b1);  // x1 = {3,2,1}
  x2 := Matrices.LU_solve(LU, pivots, b2);  // x2 = {1,0,2}

See also

Matrices.LU_solve, Matrices.solve,

Extends from Modelica.Icons.Function (Icon for a function).

Inputs

NameDescription
A[:, :]Square or rectangular matrix

Outputs

NameDescription
LU[size(A, 1), size(A, 2)]L,U factors (used with LU_solve(..))
pivots[min(size(A, 1), size(A, 2))]pivot indices (used with LU_solve(..))
infoInformation

Modelica_LinearSystems2.Math.Matrices.LU_solve Modelica_LinearSystems2.Math.Matrices.LU_solve

Solve real system of linear equations P*L*U*x=b with a b vector and an LU decomposition (from LU(..))

Information


Syntax

Matrices.LU_solve(LU, pivots, b);

Description

This function call returns the solution x of the linear systems of equations

P*L*U*x = b;

where P is a permutation matrix (implicitely defined by vector pivots), L is a lower triangular matrix with unit diagonal elements (lower trapezoidal if m > n), and U is an upper triangular matrix (upper trapezoidal if m < n). The matrices of this decomposition are computed with function Matrices.LU that returns arguments LU and pivots used as input arguments of Matrices.LU_solve. With Matrices.LU and Matrices.LU_solve it is possible to efficiently solve linear systems with different right hand side vectors. If a linear system of equations with just one right hand side vector shall be solved, it is more convenient to just use the function Matrices.solve.

If a unique solution x does not exist (since the LU decomposition is singular), an exception is raised.

The LU factorization is computed with the LAPACK function "dgetrf", i.e., by Gaussian elemination using partial pivoting with row interchanges. Vector "pivots" are the pivot indices, i.e., for 1 ≤ i ≤ min(m,n), row i of matrix A was interchanged with row pivots[i].

Example

  Real A[3,3] = [1,2,3; 
                 3,4,5;
                 2,1,4];
  Real b1[3] = {10,22,12};
  Real b2[3] = { 7,13,10};
  Real    LU[3,3];
  Integer pivots[3];
  Real    x1[3];
  Real    x2[3];
algorithm
  (LU, pivots) := Matrices.LU(A);
  x1 := Matrices.LU_solve(LU, pivots, b1);  // x1 = {3,2,1}
  x2 := Matrices.LU_solve(LU, pivots, b2);  // x2 = {1,0,2}

See also

Matrices.LU, Matrices.solve,

Extends from Modelica.Icons.Function (Icon for a function).

Inputs

NameDescription
LU[:, size(LU, 1)]L,U factors of Matrices.LU(..) for a square matrix
pivots[size(LU, 1)]Pivots indices of Matrices.LU(..)
b[size(LU, 1)]Right hand side vector of P*L*U*x=b

Outputs

NameDescription
x[size(b, 1)]Solution vector such that P*L*U*x = b

Modelica_LinearSystems2.Math.Matrices.LU_solve2

Solve real system of linear equations P*L*U*X=B with a B vector and an LU decomposition (from LU(..))

Information


Syntax

Matrices.LU_solve(LU, pivots, B);

Description

This function call returns the solution x of the linear systems of equations

P*L*U*X = B;

where P is a permutation matrix (implicitely defined by vector pivots), L is a lower triangular matrix with unit diagonal elements (lower trapezoidal if m > n), and U is an upper triangular matrix (upper trapezoidal if m < n). The matrices of this decomposition are computed with function Matrices.LU that returns arguments LU and pivots used as input arguments of Matrices.LU_solve. With Matrices.LU and Matrices.LU_solve it is possible to efficiently solve linear systems with different right hand side matrices. If a linear system of equations with just one right hand side matrix shall be solved, it is more convenient to just use the function Matrices.solve2.

If a unique solution X does not exist (since the LU decomposition is singular), an exception is raised.

The LU factorization is computed with the LAPACK function "dgetrf", i.e., by Gaussian elemination using partial pivoting with row interchanges. Vector "pivots" are the pivot indices, i.e., for 1 ≤ i ≤ min(m,n), row i of matrix A was interchanged with row pivots[i].

Example

  Real A[3,3] = [1,2,3; 
                 3,4,5;
                 2,1,4];
  Real B1[3] = [10, 20;
                22, 44;
                12, 24];
  Real B2[3] = [ 7, 14;
                13, 26;
                10, 20];
  Real    LU[3,3];
  Integer pivots[3];
  Real    X1[3,2];
  Real    X2[3,2];
algorithm
  (LU, pivots) := Matrices.LU(A);
  X1 := Matrices.LU_solve2(LU, pivots, B1);  /* X1 = [3, 6;
                                                      2, 4;
                                                      1, 2] */
  X2 := Matrices.LU_solve2(LU, pivots, B2);  /* X2 = [1, 2;
                                                      0, 0;
                                                      2, 4] */

See also

Matrices.LU, Matrices.solve2,

Inputs

NameDescription
LU[:, size(LU, 1)]L,U factors of Matrices.LU(..) for a square matrix
pivots[size(LU, 1)]Pivots indices of Matrices.LU(..)
B[size(LU, 1), :]Right hand side matrix of P*L*U*X=B

Outputs

NameDescription
X[size(B, 1), size(B, 2)]Solution matrix such that P*L*U*X = B

Modelica_LinearSystems2.Math.Matrices.lyapunov

Solution of continuous-time Lyapunov equation X*A + A'*X = C

Information


 
 
Function laypunov computes the solution X of the continuous-time Lyapunov equation
 XA + A'*X = C.
using the Schur method for Lyapunov equations proposed by Bartels and Stewart [1].

References

  [1] Bartels, R.H. and Stewart G.W.
      Algorithm 432: Solution of the matrix equation AX + XB = C.
      Comm. ACM., Vol. 15, pp. 820-826, 1972.

Inputs

NameDescription
A[:, size(A, 1)] 
C[size(A, 1), size(A, 2)] 
eps 

Outputs

NameDescription
X[size(A, 1), size(A, 2)]solution of the Riccati equation

Modelica_LinearSystems2.Math.Matrices.norm Modelica_LinearSystems2.Math.Matrices.norm

Returns the norm of a matrix

Information


Syntax

Matrices.norm(A);
Matrices.norm(A, p=2);

Description

The function call "Matrices.norm(A)" returns the 2-norm of matrix A, i.e., the largest singular value of A.
The function call "Matrices.norm(A, p)" returns the p-norm of matrix A. The only allowed values for p are

Note, for any matrices A1, A2 the following inequality holds:

Matrices.norm(A1+A2,p) ≤ Matrices.norm(A1,p) + Matrices.norm(A2,p)

Note, for any matrix A and vector v the following inequality holds:

Vectors.norm(A*v,p) ≤ Matrices.norm(A,p)*Vectors.norm(A,p)

Extends from Modelica.Icons.Function (Icon for a function).

Inputs

NameDescription
A[:, :]Input matrix
pType of p-norm (only allowed: 1, 2 or Modelica.Constants.inf)

Outputs

NameDescription
resultp-norm of matrix A

Modelica_LinearSystems2.Math.Matrices.orthogonalQ

generates a real orthogonal matrix Q which is defined as the product of IHI-ILO elementary reflectors of order N, as returned by DGEHRD

Information

   Purpose  
/*  ======= */
 
/*  DORGHR generates a real orthogonal matrix Q which is defined as the */
/*  product of IHI-ILO elementary reflectors of order N, as returned by */
/*  DGEHRD: */
 
/*  Q = H(ilo) H(ilo+1) . . . H(ihi-1). */
 
/*  Arguments */
/*  ========= */
 
/*  N       (input) INTEGER */
/*          The order of the matrix Q. N >= 0. */
 
/*  ILO     (input) INTEGER */
/*  IHI     (input) INTEGER */
/*          ILO and IHI must have the same values as in the previous call */
/*          of DGEHRD. Q is equal to the unit matrix except in the */
/*          submatrix Q(ilo+1:ihi,ilo+1:ihi). */
/*          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. */
 
/*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
/*          On entry, the vectors which define the elementary reflectors, */
/*          as returned by DGEHRD. */
/*          On exit, the N-by-N orthogonal matrix Q. */
 
/*  LDA     (input) INTEGER */
/*          The leading dimension of the array A. LDA >= max(1,N). */
 
/*  TAU     (input) DOUBLE PRECISION array, dimension (N-1) */
/*          TAU(i) must contain the scalar factor of the elementary */
/*          reflector H(i), as returned by DGEHRD. */
 
/*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK) */
/*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
 
/*  LWORK   (input) INTEGER */
/*          The dimension of the array WORK. LWORK >= IHI-ILO. */
/*          For optimum performance LWORK >= (IHI-ILO)*NB, where NB is */
/*          the optimal blocksize. */
 
/*          If LWORK = -1, then a workspace query is assumed; the routine */
/*          only calculates the optimal size of the WORK array, returns */
/*          this value as the first entry of the WORK array, and no error */
/*          message related to LWORK is issued by XERBLA. */
 
/*  INFO    (output) INTEGER */
/*          = 0:  successful exit */
/*          < 0:  if INFO = -i, the i-th argument had an illegal value */
 
/*  ===================================================================== */

Inputs

NameDescription
A[:, size(A, 1)] 
tau[size(A, 1) - 1]scalar factors of the elementary reflectors
ilolowest index where the original matrix had been Hessenbergform - ilo must have the same value as in the previous call of DGEHRD
ihihighest index where the original matrix had been Hessenbergform - ihi must have the same value as in the previous call of DGEHRD

Outputs

NameDescription
Q[size(A, 1), size(A, 2)]Orthogonal matrix as a result of elementary reflectors
info 

Modelica_LinearSystems2.Math.Matrices.printMatrix

print matrix

Inputs

NameDescription
M[:, :] 
significantDigitsNumber of significant digits that are shown
nameIndependent variable name used for printing

Outputs

NameDescription
s 

Modelica_LinearSystems2.Math.Matrices.QR

QR decomposition of a square matrix without column pivoting (A = Q*R)

Inputs

NameDescription
A[:, :]Rectangular matrix with size(A,1) >= size(A,2)

Outputs

NameDescription
Q[size(A, 1), size(A, 2)]Rectangular matrix with orthonormal columns such that Q*R=A[:,p]
R[min(size(A, 1), size(A, 2)), size(A, 2)]Square upper triangular matrix
tau[min(size(A, 1), size(A, 2))] 
Q2[:, :] 

Modelica_LinearSystems2.Math.Matrices.rsf

Computes the real Schur form (RSF) of a square matrix

Inputs

NameDescription
A[:, size(A, 1)] 

Outputs

NameDescription
S[size(A, 1), size(A, 2)] 
QZ[size(A, 1), size(A, 2)] 
alphaReal[size(A, 1)]Real part of eigenvalue=alphaReal+i*alphaImag
alphaImag[size(A, 1)]Imaginary part of eigenvalue=(alphaReal+i*alphaImag

Modelica_LinearSystems2.Math.Matrices.solve Modelica_LinearSystems2.Math.Matrices.solve

Solve real system of linear equations A*x=b with a b vector (Gaussian elemination with partial pivoting)

Information


Syntax

Matrices.solve(A,b);

Description

This function call returns the solution x of the linear system of equations

A*x = b

If a unique solution x does not exist (since A is singular), an exception is raised.

Note, the solution is computed with the LAPACK function "dgesv", i.e., by Gaussian elemination with partial pivoting.

Example

  Real A[3,3] = [1,2,3; 
                 3,4,5;
                 2,1,4];
  Real b[3] = {10,22,12};
  Real x[3];
algorithm
  x := Matrices.solve(A,b);  // x = {3,2,1}

See also

Matrices.LU, Matrices.LU_solve

Extends from Modelica.Icons.Function (Icon for a function).

Inputs

NameDescription
A[:, size(A, 1)]Matrix A of A*x = b
b[size(A, 1)]Vector b of A*x = b

Outputs

NameDescription
x[size(b, 1)]Vector x such that A*x = b

Modelica_LinearSystems2.Math.Matrices.solve2 Modelica_LinearSystems2.Math.Matrices.solve2

Solve real system of linear equations A*X=B with a B matrix (Gaussian elemination with partial pivoting)

Information


Syntax

Matrices.solve2(A,b);

Description

This function call returns the solution X of the linear system of equations

A*X = B

If a unique solution X does not exist (since A is singular), an exception is raised.

Note, the solution is computed with the LAPACK function "dgesv", i.e., by Gaussian elemination with partial pivoting.

Example

  Real A[3,3] = [1,2,3; 
                 3,4,5;
                 2,1,4];
  Real B[3,2] = [10, 20;
                 22, 44;
                 12, 24];
  Real X[3,2];
algorithm
  (LU, pivots) := Matrices.LU(A);
  X := Matrices.solve2(A, B1);  /* X = [3, 6;
                                        2, 4;
                                        1, 2] */

See also

Matrices.LU, Matrices.LU_solve2

Extends from Modelica.Icons.Function (Icon for a function).

Inputs

NameDescription
A[:, size(A, 1)]Matrix A of A*X = B
B[size(A, 1), :]Matrix B of A*X = B

Outputs

NameDescription
X[size(B, 1), size(B, 2)]Matrix X such that A*X = B

Modelica_LinearSystems2.Math.Matrices.toUpperHessenberg

transform a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation: Q' * A * Q = H

Information


This function uses DGEHRD-LAPACK routine to calculate a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation:  Q' * A * Q = H.

Inputs are

Outputs are

 
 
   =====================================================================  
    -- LAPACK routine (version 3.0) --  
      Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,  
      Courant Institute, Argonne National Lab, and Rice University  
      June 30, 1999  
 
  
   Purpose  
   =======  
 
   DGEHRD reduces a real general matrix A to upper Hessenberg form H by  
   an orthogonal similarity transformation:  Q' * A * Q = H .  
 
   Arguments  
   =========  
 
   N       (input) INTEGER  
           The order of the matrix A.  N >= 0.  
 
   ILO     (input) INTEGER  
   IHI     (input) INTEGER  
           It is assumed that A is already upper triangular in rows  
           and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally  
           set by a previous call to DGEBAL; otherwise they should be  
           set to 1 and N respectively. See Further Details.  
           1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.  
 
   A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)  
           On entry, the N-by-N general matrix to be reduced.  
           On exit, the upper triangle and the first subdiagonal of A  
           are overwritten with the upper Hessenberg matrix H, and the  
           elements below the first subdiagonal, with the array TAU,  
           represent the orthogonal matrix Q as a product of elementary  
           reflectors. See Further Details.  
 
   LDA     (input) INTEGER  
           The leading dimension of the array A.  LDA >= max(1,N).  
 
   TAU     (output) DOUBLE PRECISION array, dimension (N-1)  
           The scalar factors of the elementary reflectors (see Further  
           Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to  
           zero.  
 
   WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)  
           On exit, if INFO = 0, WORK(1) returns the optimal LWORK.  
 
   LWORK   (input) INTEGER  
           The length of the array WORK.  LWORK >= max(1,N).  
           For optimum performance LWORK >= N*NB, where NB is the  
           optimal blocksize.  
 
           If LWORK = -1, then a workspace query is assumed; the routine  
           only calculates the optimal size of the WORK array, returns  
           this value as the first entry of the WORK array, and no error  
           message related to LWORK is issued by XERBLA.  
 
   INFO    (output) INTEGER  
           = 0:  successful exit  
           < 0:  if INFO = -i, the i-th argument had an illegal value.  
 
   Further Details  
   ===============  
 
   The matrix Q is represented as a product of (ihi-ilo) elementary  
   reflectors  
 
      Q = H(ilo) H(ilo+1) . . . H(ihi-1).  
 
   Each H(i) has the form  
 
      H(i) = I - tau * v * v'  
 
   where tau is a real scalar, and v is a real vector with  
   v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on  
   exit in A(i+2:ihi,i), and tau in TAU(i).  
 
   The contents of A are illustrated by the following example, with  
   n = 7, ilo = 2 and ihi = 6:  
 
   on entry,                        on exit,  
 
   ( a   a   a   a   a   a   a )    (  a   a   h   h   h   h   a )  
   (     a   a   a   a   a   a )    (      a   h   h   h   h   a )  
   (     a   a   a   a   a   a )    (      h   h   h   h   h   h )  
   (     a   a   a   a   a   a )    (      v2  h   h   h   h   h )  
   (     a   a   a   a   a   a )    (      v2  v3  h   h   h   h )  
   (     a   a   a   a   a   a )    (      v2  v3  v4  h   h   h )  
   (                         a )    (                          a )  
 
   where a denotes an element of the original matrix A, h denotes a  
   modified element of the upper Hessenberg matrix H, and vi denotes an  
   element of the vector defining H(i).  
 
   =====================================================================  

Inputs

NameDescription
A[:, size(A, 1)] 
ilolowest index where the original matrix had been Hessenbergform
ihi 

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[size(A, 1) - 1] 
info 

Modelica_LinearSystems2.Math.Matrices.trace Modelica_LinearSystems2.Math.Matrices.trace

tarce(A) is the sum of the diagonal elements of A

Information

Extends from Modelica.Icons.Function (Icon for a function).

Inputs

NameDescription
A[:, size(A, 1)] 

Outputs

NameDescription
result 

HTML-documentation generated by Dymola Tue Sep 08 18:52:57 2009.