Modelica.Math.Matrices

Information

Library content

This library provides functions operating on matrices:

Most functions are solely an interface to the external LAPACK library (http://www.netlib.org/lapack). The details of this library are described in:

Anderson E., Bai Z., Bischof C., Blackford S., Demmel J., Dongarra J., Du Croz J., Greenbaum A., Hammarling S., McKenney A., and Sorensen D.:

Lapack Users' Guide. Third Edition, SIAM, 1999.

See also

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

Package Content

Name	Description
isEqual	Compare whether two Real matrices are identical
norm	Returns the norm of a matrix
sort	Sort rows or columns of matrix in ascending or descending order
solve	Solve real system of linear equations A*x=b with a b vector (Gaussian elemination with partial pivoting)
solve2	Solve real system of linear equations A*X=B with a B matrix (Gaussian elemination with partial pivoting)
leastSquares	Solve overdetermined or underdetermined real system of linear equations A*x=b in a least squares sense (A may be rank deficient)
equalityLeastSquares	Solve a linear equality constrained least squares problem
LU	LU decomposition of square or rectangular matrix
LU_solve	Solve real system of linear equations P*L*U*x=b with a b vector and an LU decomposition (from LU(..))
LU_solve2	Solve real system of linear equations P*L*U*X=B with a B matrix and an LU decomposition (from LU(..))
QR	QR decomposition of a square matrix with column pivoting (A(:,p) = Q*R)
eigenValues	Compute eigenvalues and eigenvectors for a real, nonsymmetric matrix
eigenValueMatrix	Return real valued block diagonal matrix J of eigenvalues of matrix A (A=V*J*Vinv)
singularValues	Compute singular values and left and right singular vectors
det	Determinant of a matrix (computed by LU decomposition)
inv	Inverse of a matrix (try to avoid, use function solve(..) instead)
rank	Rank of a matrix (computed with singular values)
balance	Balancing of matrix A to improve the condition of A
exp	Compute the exponential of a matrix by adaptive Taylor series expansion with scaling and balancing
integralExp	Computation of the transition-matrix phi and its integral gamma
integralExpT	Computation of the transition-matrix phi and the integral gamma and gamma1
LAPACK	Interface to LAPACK library (should usually not directly be used but only indirectly via Modelica.Math.Matrices)

Modelica.Math.Matrices.isEqual

Information

Syntax

Matrices.isEqual(M1, M2);
Matrices.isEqual(M1, M2, eps=0);

Description

The function call "Matrices.isEqual(M1, M2)" returns true, if the two Real matrices M1 and M2 have the same dimensions and the same elements. Otherwise the function returns false. Two elements e1 and e2 of the two matrices are checked on equality by the test "abs(e1-e2) ≤ eps", where "eps" can be provided as third argument of the function. Default is "eps = 0".

Example

  Real A1[2,2] = [1,2; 3,4];
  Real A2[3,2] = [1,2; 3,4; 5,6];
  Real A3[2,2] = [1,2, 3,4.0001];
  Boolean result;
algorithm
  result := Matrices.isEqual(M1,M2);     // = false
  result := Matrices.isEqual(M1,M3);     // = false
  result := Matrices.isEqual(M1,M1);     // = true
  result := Matrices.isEqual(M1,M3,0.1); // = true

See also

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

Inputs

Type	Name	Default	Description
Real	M1[:, :]		First matrix
Real	M2[:, :]		Second matrix (may have different size as M1
Real	eps	0	Two elements e1 and e2 of the two matrices are identical if abs(e1-e2) <= eps

Outputs

Type	Name	Description
Boolean	result	= true, if matrices have the same size and the same elements

Modelica definition

function isEqual "Compare whether two Real matrices are identical"
  extends Modelica.Icons.Function;
  input Real M1[:, :] "First matrix";
  input Real M2[:, :] "Second matrix (may have different size as M1";
  input Real eps(min=0) = 0 
    "Two elements e1 and e2 of the two matrices are identical if abs(e1-e2) <= eps";
  output Boolean result 
    "= true, if matrices have the same size and the same elements";

protected 
  Integer nrow=size(M1, 1) "Number of rows of matrix M1";
  Integer ncol=size(M1, 2) "Number of columns of matrix M1";
  Integer i=1;
  Integer j;
algorithm 
  result := false;
  if size(M2, 1) == nrow and size(M2, 2) == ncol then
    result := true;
    while i <= nrow loop
      j := 1;
      while j <= ncol loop
        if abs(M1[i, j] - M2[i, j]) > eps then
          result := false;
          i := nrow;
          j := ncol;
        end if;
        j := j + 1;
      end while;
      i := i + 1;
    end while;
  end if;

end isEqual;

Modelica.Math.Matrices.norm

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

• "p=1": the largest column sum of A
• "p=2": the largest singular value of A
• "p=Modelica.Constants.inf": the largest row sum of A

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

Type	Name	Default	Description
Real	A[:, :]		Input matrix
Real	p	2	Type of p-norm (only allowed: 1, 2 or Modelica.Constants.inf)

Outputs

Type	Name	Description
Real	result	p-norm of matrix A

Modelica definition

function norm "Returns the norm of a matrix"
  extends Modelica.Icons.Function;
  input Real A[:, :] "Input matrix";
  input Real p(min=1) = 2 
    "Type of p-norm (only allowed: 1, 2 or Modelica.Constants.inf)";
  output Real result=0.0 "p-norm of matrix A";

algorithm 
  if p == 1 then
    // column sum norm
    for i in 1:size(A, 2) loop
      result := max(result, sum(abs(A[:, i])));
    end for;
  elseif p == 2 then
    // largest singular value
    result := max(singularValues(A));
  elseif p == Modelica.Constants.inf then
    // row sum norm
    for i in 1:size(A, 1) loop
      result := max(result, sum(abs(A[i, :])));
    end for;
  else
    assert(false, "Optional argument \"p\" of function \"norm\" must be
1, 2 or Modelica.Constants.inf");
  end if;
end norm;

Modelica.Math.Matrices.sort

Information

Syntax

           sorted_M = Matrices.sort(M);
(sorted_M, indices) = Matrices.sort(M, sortRows=true, ascending=true);

Description

Function sort(..) sorts the rows of a Real matrix M in ascending order and returns the result in sorted_M. If the optional argument "sortRows" is false, the columns of the matrix are sorted. If the optional argument "ascending" is false, the rows or columns are sorted in descending order. In the optional second output argument, the indices of the sorted rows or columns with respect to the original matrix are given, such that

   sorted_M = if sortedRow then M[indices,:] else M[:,indices];

Example

  (M2, i2) := Matrices.sort([2, 1,  0;
                             2, 0, -1]);
       -> M2 = [2, 0, -1;
                2, 1, 0 ];
          i2 = {2,1};

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

Inputs

Type	Name	Default	Description
Real	M[:, :]		Matrix to be sorted
Boolean	sortRows	true	= true if rows are sorted, otherwise columns
Boolean	ascending	true	= true if ascending order, otherwise descending order

Outputs

Type	Name	Description
Real	sorted_M[size(M, 1), size(M, 2)]	Sorted matrix
Integer	indices[if sortRows then size(M, 1) else size(M, 2)]	sorted_M = if sortRows then M[indices,:] else M[:,indices]

Modelica definition

function sort 
  "Sort rows or columns of matrix in ascending or descending order"
  extends Modelica.Icons.Function;
  input Real M[:,:] "Matrix to be sorted";
  input Boolean sortRows = true "= true if rows are sorted, otherwise columns";
  input Boolean ascending = true 
    "= true if ascending order, otherwise descending order";
  output Real sorted_M[size(M,1), size(M,2)] = M "Sorted matrix";
  output Integer indices[if sortRows then size(M,1) else size(M,2)] 
    "sorted_M = if sortRows then M[indices,:] else M[:,indices]";

  /* shellsort algorithm; should be improved later */
protected 
  Integer gap;
  Integer i;
  Integer j;
  Real wM2[size(M,2)];
  Integer wi;
  Integer nM1 = size(M,1);
  Boolean swap;
  Real sorted_MT[size(M,2), size(M,1)];

encapsulated function greater "Compare whether vector v1 > v2"
    import Modelica;
  extends Modelica.Icons.Function;
    import Modelica.Utilities.Types.Compare;
  input Real v1[:];
  input Real v2[size(v1,1)];
  output Boolean result;
  protected 
  Integer n = size(v1,1);
  Integer i=1;
algorithm 
  result := false;
  while i <= n loop
     if v1[i] > v2[i] then
        result := true;
        i := n;
     elseif v1[i] < v2[i] then
        i := n;
     end if;
     i := i+1;
  end while;
end greater;

encapsulated function less "Compare whether vector v1 < v2"
    import Modelica;
  extends Modelica.Icons.Function;
    import Modelica.Utilities.Types.Compare;
  input Real v1[:];
  input Real v2[size(v1,1)];
  output Boolean result;
  protected 
  Integer n = size(v1,1);
  Integer i=1;
algorithm 
  result := false;
  while i <= n loop
     if v1[i] < v2[i] then
        result := true;
        i := n;
     elseif v1[i] > v2[i] then
        i := n;
     end if;
     i := i+1;
  end while;
end less;
algorithm 
  if not sortRows then
      (sorted_MT,indices) := sort(transpose(M), ascending=ascending);
     sorted_M :=transpose(sorted_MT);
  else
     indices :=1:size(M, 1);
     gap := div(nM1,2);
     while gap > 0 loop
        i := gap;
        while i < nM1 loop
           j := i-gap;
           if j>=0 then
              if ascending then
                 swap := greater(sorted_M[j+1,:], sorted_M[j+gap+1,:]);
              else
                 swap := less(sorted_M[j+1,:], sorted_M[j+gap+1,:]);
              end if;
           else
              swap := false;
           end if;

           while swap loop
              wM2 := sorted_M[j+1,:];
              wi := indices[j+1];
              sorted_M[j+1,:] := sorted_M[j+gap+1,:];
              sorted_M[j+gap+1,:] := wM2;
              indices[j+1] := indices[j+gap+1];
              indices[j+gap+1] := wi;
              j := j - gap;
              if j >= 0 then
                 if ascending then
                    swap := greater(sorted_M[j+1,:], sorted_M[j+gap+1,:]);
                 else
                    swap := less(sorted_M[j+1,:], sorted_M[j+gap+1,:]);
                 end if;
              else
                 swap := false;
              end if;
           end while;
           i := i + 1;
        end while;
        gap := div(gap,2);
     end while;
  end if;
end sort;

Modelica.Math.Matrices.solve

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

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

Inputs

Type	Name	Default	Description
Real	A[:, size(A, 1)]		Matrix A of A*x = b
Real	b[size(A, 1)]		Vector b of A*x = b

Outputs

Type	Name	Description
Real	x[size(b, 1)]	Vector x such that A*x = b

Modelica definition

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

  extends Modelica.Icons.Function;
  input Real A[:, size(A, 1)] "Matrix A of A*x = b";
  input Real b[size(A, 1)] "Vector b of A*x = b";
  output Real x[size(b, 1)] "Vector x such that A*x = b";

protected 
  Integer info;
algorithm 
  (x,info) := LAPACK.dgesv_vec(A, b);
  assert(info == 0, "Solving a linear system of equations with function
\"Matrices.solve\" is not possible, because the system has either
no or infinitely many solutions (A is singular).");
end solve;

Modelica.Math.Matrices.solve2

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

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

Inputs

Type	Name	Default	Description
Real	A[:, size(A, 1)]		Matrix A of A*X = B
Real	B[size(A, 1), :]		Matrix B of A*X = B

Outputs

Type	Name	Description
Real	X[size(B, 1), size(B, 2)]	Matrix X such that A*X = B

Modelica definition

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

  extends Modelica.Icons.Function;
  input Real A[:, size(A, 1)] "Matrix A of A*X = B";
  input Real B[size(A, 1),:] "Matrix B of A*X = B";
  output Real X[size(B, 1), size(B,2)] "Matrix X such that A*X = B";

protected 
  Integer info;
algorithm 
  (X,info) := LAPACK.dgesv(A, B);
  assert(info == 0, "Solving a linear system of equations with function
\"Matrices.solve2\" is not possible, because the system has either
no or infinitely many solutions (A is singular).");
end solve2;

Modelica.Math.Matrices.leastSquares

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

Type	Name	Default	Description
Real	A[:, :]		Matrix A
Real	b[size(A, 1)]		Vector b

Outputs

Type	Name	Description
Real	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 definition

function leastSquares 
  "Solve overdetermined or underdetermined real system of linear equations A*x=b in a least squares sense (A may be rank deficient)"
  extends Modelica.Icons.Function;
  input Real A[:, :] "Matrix A";
  input Real b[size(A, 1)] "Vector b";
  output Real 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)";

protected 
  Integer info;
  Integer rank;
  Real xx[max(size(A,1),size(A,2))];
algorithm 
  (xx,info,rank) := LAPACK.dgelsx_vec(A, b, 100*Modelica.Constants.eps);
  x := xx[1:size(A,2)];
  assert(info == 0, "Solving an overdetermined or underdetermined linear system of
equations with function \"Matrices.leastSquares\" failed.");
end leastSquares;

Modelica.Math.Matrices.equalityLeastSquares

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

Type	Name	Default	Description
Real	A[:, :]		Minimize |A*x - a|^2
Real	a[size(A, 1)]
Real	B[:, size(A, 2)]		subject to B*x=b
Real	b[size(B, 1)]

Outputs

Type	Name	Description
Real	x[size(A, 2)]	solution vector

Modelica definition

function equalityLeastSquares 
  "Solve a linear equality constrained least squares problem"
  extends Modelica.Icons.Function;
  input Real A[:,:] "Minimize |A*x - a|^2";
  input Real a[size(A,1)];
  input Real B[:,size(A,2)] "subject to B*x=b";
  input Real b[size(B,1)];
  output Real x[size(A,2)] "solution vector";

protected 
  Integer info;
algorithm 
  assert(size(A,2) >= size(B,1) and size(A,2) <= size(A,1) + size(B,1),
         "It is required that size(B,1) <= size(A,2) <= size(A,1) + size(B,1)\n" +
         "This relationship is not fulfilled, since the matrices are declared as:\n" +
         "  A[" + String(size(A,1)) + "," + String(size(A,2)) + "], B[" +
         String(size(B,1)) + "," + String(size(B,2)) + "]\n");

  (x, info) := LAPACK.dgglse_vec(A, a, B, b);

  assert(info == 0, "Solving a linear equality-constrained least squares problem
with function \"Matrices.equalityLeastSquares\" failed.");
end equalityLeastSquares;

Modelica.Math.Matrices.LU

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

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

Inputs

Type	Name	Default	Description
Real	A[:, :]		Square or rectangular matrix

Outputs

Type	Name	Description
Real	LU[size(A, 1), size(A, 2)]	L,U factors (used with LU_solve(..))
Integer	pivots[min(size(A, 1), size(A, 2))]	pivot indices (used with LU_solve(..))
Integer	info	Information

Modelica definition

function LU "LU decomposition of square or rectangular matrix"
  extends Modelica.Icons.Function;
  input Real A[:, :] "Square or rectangular matrix";
  output Real LU[size(A, 1), size(A,2)] = A 
    "L,U factors (used with LU_solve(..))";
  output Integer pivots[min(size(A, 1), size(A,2))] 
    "pivot indices (used with LU_solve(..))";
  output Integer info "Information";
  external "FORTRAN 77" dgetrf(size(A, 1), size(A, 2), LU, size(A, 1), pivots, info);

end LU;

Modelica.Math.Matrices.LU_solve

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

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

Inputs

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

Outputs

Type	Name	Description
Real	x[size(b, 1)]	Solution vector such that P*L*U*x = b

Modelica definition

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

  extends Modelica.Icons.Function;
  input Real LU[:, size(LU,1)] 
    "L,U factors of Matrices.LU(..) for a square matrix";
  input Integer pivots[size(LU, 1)] "Pivots indices of Matrices.LU(..)";
  input Real b[size(LU, 1)] "Right hand side vector of P*L*U*x=b";
  output Real x[size(b, 1)] "Solution vector such that P*L*U*x = b";

algorithm 
  for i in 1:size(LU,1) loop
       assert(LU[i,i] <> 0, "Solving a linear system of equations with function
\"Matrices.LU_solve\" is not possible, since the LU decomposition
is singular, i.e., no unique solution exists.");
  end for;
  x := LAPACK.dgetrs_vec(LU, pivots, b);
end LU_solve;

Modelica.Math.Matrices.LU_solve2

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_solve2. With Matrices.LU and Matrices.LU_solve2 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

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

Inputs

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

Outputs

Type	Name	Description
Real	X[size(B, 1), size(B, 2)]	Solution matrix such that P*L*U*X = B

Modelica definition

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

  extends Modelica.Icons.Function;
  input Real LU[:, size(LU,1)] 
    "L,U factors of Matrices.LU(..) for a square matrix";
  input Integer pivots[size(LU, 1)] "Pivots indices of Matrices.LU(..)";
  input Real B[size(LU, 1),:] "Right hand side matrix of P*L*U*X=B";
  output Real X[size(B, 1), size(B,2)] "Solution matrix such that P*L*U*X = B";

algorithm 
  for i in 1:size(LU,1) loop
       assert(LU[i,i] <> 0, "Solving a linear system of equations with function
\"Matrices.LU_solve\" is not possible, since the LU decomposition
is singular, i.e., no unique solution exists.");
  end for;
  X := Modelica.Math.Matrices.LAPACK.dgetrs(
    LU,
    pivots,
    B);
end LU_solve2;

Modelica.Math.Matrices.QR

Information

Syntax

(Q,R,p) = Matrices.QR(A);

Description

This function returns the QR decomposition of a rectangular matrix A (the number of columns of A must be less than or equal to the number of rows):

Q*R = A[:,p]

where Q is a rectangular matrix that has orthonormal columns and has the same size as A (Q^TQ=I), R is a square, upper triangular matrix and p is a permutation vector. Matrix R has the following important properties:

• The absolute value of a diagonal element of R is the largest value in this row, i.e., abs(R[i,i]) ≥ abs(R[i,j]).
• The diagonal elements of R are sorted according to size, such that the largest absolute value is abs(R[1,1]) and abs(R[i,i]) ≥ abs(R[j,j]) with i < j.

This means that if abs(R[i,i]) ≤ ε then abs(R[j,k]) ≤ ε for j ≥ i, i.e., the i-th row up to the last row of R have small elements and can be treated as being zero. This allows to, e.g., estimate the row-rank of R (which is the same row-rank as A). Furthermore, R can be partitioned in two parts

   A[:,p] = Q * [R1, R2;
                 0,  0]

where R₁ is a regular, upper triangular matrix.

Note, the solution is computed with the LAPACK functions "dgeqpf" and "dorgqr", i.e., by Housholder transformations with column pivoting. If Q is not needed, the function may be called as: (,R,p) = QR(A).

Example

  Real A[3,3] = [1,2,3;
                 3,4,5;
                 2,1,4];
  Real R[3,3];
algorithm
  (,R) := Matrices.QR(A);  // R = [-7.07.., -4.24.., -3.67..;
                                    0     , -1.73.., -0.23..;
                                    0     ,  0     ,  0.65..];

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

Inputs

Type	Name	Default	Description
Real	A[:, :]		Rectangular matrix with size(A,1) >= size(A,2)

Outputs

Type	Name	Description
Real	Q[size(A, 1), size(A, 2)]	Rectangular matrix with orthonormal columns such that Q*R=A[:,p]
Real	R[size(A, 2), size(A, 2)]	Square upper triangular matrix
Integer	p[size(A, 2)]	Column permutation vector

Modelica definition

function QR 
  "QR decomposition of a square matrix with column pivoting (A(:,p) = Q*R)"

  extends Modelica.Icons.Function;
  input Real A[:, :] "Rectangular matrix with size(A,1) >= size(A,2)";
  output Real Q[size(A, 1), size(A, 2)] 
    "Rectangular matrix with orthonormal columns such that Q*R=A[:,p]";
  output Real R[size(A, 2), size(A, 2)] "Square upper triangular matrix";
  output Integer p[size(A, 2)] "Column permutation vector";

protected 
  Integer nrow=size(A, 1);
  Integer ncol=size(A, 2);
  Real tau[ncol];
algorithm 
  assert(nrow >= ncol, "\nInput matrix A[" + String(nrow) + "," + String(ncol) + "] has more columns as rows.
This is not allowed when calling Modelica.Matrices.QR(A).");
  (Q,tau,p) := LAPACK.dgeqpf(A);

  // determine R
  R := zeros(ncol,ncol);
  for i in 1:ncol loop
    for j in i:ncol loop
      R[i, j] := Q[i,j];
    end for;
  end for;

  // if isPresent(Q) then (not yet supported by Dymola)
  Q := LAPACK.dorgqr(Q, tau);
end QR;

Modelica.Math.Matrices.eigenValues

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; in some cases, an inversion is also possible if some eigenvalues are the same).

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

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

Inputs

Type	Name	Default	Description
Real	A[:, size(A, 1)]		Matrix

Outputs

Type	Name	Description
Real	eigenvalues[size(A, 1), 2]	Eigenvalues of matrix A (Re: first column, Im: second column)
Real	eigenvectors[size(A, 1), size(A, 2)]	Real-valued eigenvector matrix

Modelica definition

function eigenValues 
  "Compute eigenvalues and eigenvectors for a real, nonsymmetric matrix"

  extends Modelica.Icons.Function;
  input Real A[:, size(A, 1)] "Matrix";
  output Real eigenvalues[size(A, 1), 2] 
    "Eigenvalues of matrix A (Re: first column, Im: second column)";
  output Real eigenvectors[size(A,1), size(A,2)] 
    "Real-valued eigenvector matrix";

protected 
  Integer info;
  // replace with "isPresent(..)" if supported by Dymola
  Boolean onlyEigenvalues = false;
algorithm 
if size(A,1) > 0 then
  if onlyEigenvalues then
     (eigenvalues[:, 1],eigenvalues[:, 2],info) := LAPACK.dgeev_eigenValues(A);
     eigenvectors :=zeros(size(A, 1), size(A, 1));
  else
     (eigenvalues[:, 1],eigenvalues[:, 2],eigenvectors, info) := LAPACK.dgeev(A);
  end if;
  assert(info == 0, "Calculating the eigen values with function
\"Matrices.eigenvalues\" is not possible, since the
numerical algorithm does not converge.");
end if;
end eigenValues;

Modelica.Math.Matrices.eigenValueMatrix

Information

Syntax

Matrices.eigenValueMatrix(eigenvalues);

Description

The function call returns a block diagonal matrix J from the the two-column matrix eigenvalues (computed by function Matrices.eigenValues). Matrix eigenvalues must have the real part of the eigenvalues in the first column and the imaginary part in the second column. If an eigenvalue i has a vanishing imaginary part, then J[i,i] = eigenvalues[i,1], i.e., the diagonal element of J is the real eigenvalue. Otherwise, eigenvalue i and conjugate complex eigenvalue i+1 are used to construct a 2 by 2 diagonal block of J:

  J[i  , i]   := eigenvalues[i,1];
  J[i  , i+1] := eigenvalues[i,2];
  J[i+1, i]   := eigenvalues[i+1,2];
  J[i+1, i+1] := eigenvalues[i+1,1];

See also

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

Inputs

Type	Name	Default	Description
Real	eigenValues[:, 2]		Eigen values from function eigenValues(..) (Re: first column, Im: second column)

Outputs

Type	Name	Description
Real	J[size(eigenValues, 1), size(eigenValues, 1)]	Real valued block diagonal matrix with eigen values (Re: 1x1 block, Im: 2x2 block)

Modelica definition

function eigenValueMatrix 
  "Return real valued block diagonal matrix J of eigenvalues of matrix A (A=V*J*Vinv)"

  extends Modelica.Icons.Function;
  input Real eigenValues[:, 2] 
    "Eigen values from function eigenValues(..) (Re: first column, Im: second column)";
  output Real J[size(eigenValues, 1), size(eigenValues, 1)] 
    "Real valued block diagonal matrix with eigen values (Re: 1x1 block, Im: 2x2 block)";

protected 
  Integer n=size(eigenValues, 1);
  Integer i;
algorithm 
  J := zeros(n, n);
  i := 1;
  while i <= n loop
    if eigenValues[i, 2] == 0 then
      J[i, i] := eigenValues[i, 1];
      i := i + 1;
    else
      J[i, i] := eigenValues[i, 1];
      J[i, i + 1] := eigenValues[i, 2];
      J[i + 1, i] := eigenValues[i + 1, 2];
      J[i + 1, i + 1] := eigenValues[i + 1, 1];
      i := i + 2;
    end if;
  end while;
end eigenValueMatrix;

Modelica.Math.Matrices.singularValues

Information

Syntax

         sigma = Matrices.singularValues(A);
(sigma, U, VT) = Matrices.singularValues(A);

Description

This function computes the singular values and optionally the singular vectors of matrix A. Basically the singular value decomposition of A is computed, i.e.,

A = U S VT
  = U*Sigma*VT

where U and V are orthogonal matrices (UU^T=I, VV^T=I). S = diag(s_i) has the same size as matrix A with nonnegative diagonal elements in decreasing order and with all other elements zero (s₁ is the largest element). The function returns the singular values s_i in vector sigma and the orthogonal matrices in matrices U and V.

Example

  A = [1, 2,  3,  4;
       3, 4,  5, -2;
      -1, 2, -3,  5];
  (sigma, U, VT) = singularValues(A);
  results in:
     sigma = {8.33, 6.94, 2.31};
  i.e.
     Sigma = [8.33,    0,    0, 0;
                 0, 6.94,    0, 0;
                 0,    0, 2.31, 0]

See also

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

Inputs

Type	Name	Default	Description
Real	A[:, :]		Matrix

Outputs

Type	Name	Description
Real	sigma[min(size(A, 1), size(A, 2))]	Singular values
Real	U[size(A, 1), size(A, 1)]	Left orthogonal matrix
Real	VT[size(A, 2), size(A, 2)]	Transposed right orthogonal matrix

Modelica definition

function singularValues 
  "Compute singular values and left and right singular vectors"
  extends Modelica.Icons.Function;
  input Real A[:, :] "Matrix";
  output Real sigma[min(size(A, 1), size(A, 2))] "Singular values";
  output Real U[size(A, 1), size(A, 1)]=identity(size(A, 1)) 
    "Left orthogonal matrix";
  output Real VT[size(A, 2), size(A, 2)]=identity(size(A, 2)) 
    "Transposed right orthogonal matrix ";

protected 
  Integer info;
  Integer n=min(size(A, 1), size(A, 2)) "Number of singular values";
algorithm 
if n>0 then
  (sigma,U,VT,info) := Modelica.Math.Matrices.LAPACK.dgesvd(A);
  assert(info == 0, "The numerical algorithm to compute the
singular value decomposition did not converge");
end if;
end singularValues;

Modelica.Math.Matrices.det

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

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

Inputs

Type	Name	Default	Description
Real	A[:, size(A, 1)]

Outputs

Type	Name	Description
Real	result	Determinant of matrix A

Modelica definition

function det "Determinant of a matrix (computed by LU decomposition)"

  extends Modelica.Icons.Function;
  input Real A[:, size(A, 1)];
  output Real result "Determinant of matrix A";
protected 
  Real LU[size(A,1),size(A,1)];
  Integer pivots[size(A,1)];

algorithm 
  (LU,pivots) := Matrices.LU(A);
  result:=product(LU[i,i] for i in 1:size(A,1))*
    product(if pivots[i]==i then 1 else -1 for i in 1:size(pivots,1));
end det;

Modelica.Math.Matrices.inv

Information

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

Inputs

Type	Name	Default	Description
Real	A[:, size(A, 1)]

Outputs

Type	Name	Description
Real	invA[size(A, 1), size(A, 2)]	Inverse of matrix A

Modelica definition

function inv 
  "Inverse of a matrix (try to avoid, use function solve(..) instead)"
  extends Modelica.Icons.Function;
  input Real A[:, size(A, 1)];
  output Real invA[size(A, 1), size(A, 2)] "Inverse of matrix A";
protected 
  Integer info;
  Integer pivots[size(A, 1)] "Pivot vector";
  Real LU[size(A, 1), size(A, 2)] "LU factors of A";
algorithm 
  (LU,pivots,info) := LAPACK.dgetrf(A);

  assert(info == 0, "Calculating an inverse matrix with function
\"Matrices.inv\" is not possible, since matrix A is singular.");

  invA := LAPACK.dgetri(LU, pivots);

end inv;

Modelica.Math.Matrices.rank

Information

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

Inputs

Type	Name	Default	Description
Real	A[:, :]		Matrix
Real	eps	0	If eps > 0, the singular values are checked against eps; otherwise eps=max(size(A))*norm(A)*Modelica.Constants.eps is used

Outputs

Type	Name	Description
Integer	result	Rank of matrix A

Modelica definition

function rank "Rank of a matrix (computed with singular values)"
  extends Modelica.Icons.Function;
  input Real A[:, :] "Matrix";
  input Real eps=0 
    "If eps > 0, the singular values are checked against eps; otherwise eps=max(size(A))*norm(A)*Modelica.Constants.eps is used";
  output Integer result "Rank of matrix A";

protected 
  Integer n=min(size(A, 1), size(A, 2));
  Integer i=n;
  Real sigma[n];
  Real eps2;
algorithm 
  result := 0;
  if n > 0 then
    sigma := Modelica.Math.Matrices.singularValues(A);
    eps2 := if eps > 0 then eps else max(size(A))*sigma[1]*Modelica.Constants.eps;
    while i > 0 loop
      if sigma[i] > eps2 then
        result := i;
        i := 0;
      end if;
      i := i - 1;
    end while;
  end if;
end rank;

Modelica.Math.Matrices.balance

Information

The function transformates the matrix A, so that the norm of the i-th column is nearby the i-th row. (D,B)=Matrices.balance(A) returns a vector D, such that B=inv(diagonal(D))*A*diagonal(D) has better condition. The elements of D are multiples of 2. Balancing attempts to make the norm of each row equal to the norm of the belonging column.
Balancing is used to minimize roundoff errors inducted through large matrix calculations like Taylor-series approximation or computation of eigenvalues.

       - A = [1, 10,  1000; .01,  0,  10; .005,  .01,  10]
       - Matrices.norm(A, 1);
         = 1020.0
       - (T,B)=Matrices.balance(A)
       - T
         = {256, 16, 0.5}
       - B
         =  [1,     0.625,   1.953125;
             0.16,  0,       0.3125;
             2.56,  0.32,   10.0]
       - Matrices.norm(B, 1);
         = 12.265625

The Algorithm is taken from

H. D. Joos, G. Grbel:

RASP'91 Regulator Analysis and Synthesis Programs
DLR - Control Systems Group 1991

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

Inputs

Type	Name	Default	Description
Real	A[:, size(A, 1)]

Outputs

Type	Name	Description
Real	D[size(A, 1)]	diagonal(D)=T is transformation matrix, such that T*A*inv(T) has smaller condition as A
Real	B[size(A, 1), size(A, 1)]	Balanced matrix (= diagonal(D)*A*inv(diagonal(D)))

Modelica definition

function balance 
  "Balancing of matrix A to improve the condition of A"
  extends Modelica.Icons.Function;
  input Real A[:, size(A, 1)];
  output Real D[size(A, 1)] "diagonal(D)=T is transformation matrix, such that
          T*A*inv(T) has smaller condition as A";
  output Real B[size(A, 1), size(A, 1)] 
    "Balanced matrix (= diagonal(D)*A*inv(diagonal(D)))";
protected 
  Integer na=size(A, 1);
  Integer radix=2 "Radix of exponent representation must be 'radix'
          or a multiple of 'radix'";
  Integer radix2=radix*radix;
  Boolean noconv=true;
  Integer i=1;
  Integer j=1;
  Real CO;
  Real RO;
  Real G;
  Real F;
  Real S;
  /*auxiliary variables*/

algorithm 
  // B = inv(D)*A*D, so that cond(B)<=cond(A)
  D := ones(na);
  B := A;
  while noconv loop
    noconv := false;
    for i in 1:na loop
      CO := sum(abs(B[:, i])) - abs(B[i, i]);
      RO := sum(abs(B[i, :])) - abs(B[i, i]);
      G := RO/radix;
      F := 1;
      S := CO + RO;
      while not (CO >= G or CO == 0) loop
        F := F*radix;
        CO := CO*radix2;
      end while;
      G := RO*radix;
      while not (CO < G or RO == 0) loop
        F := F/radix;
        CO := CO/radix2;
      end while;
      if not ((CO + RO)/F >= 0.95*S) then
        G := 1/F;
        D[i] := D[i]*F;
        B[i, :] := B[i, :]*G;
        B[:, i] := B[:, i]*F;
        noconv := true;
      end if;
    end for;
  end while;
end balance;

Modelica.Math.Matrices.exp

Information

This function computes

                            (AT)^2   (AT)^3
      Φ = e^(AT) = I + AT + ------ + ------ + ....
                              2!       3!

where e=2.71828..., A is an n x n matrix with real elements and T is a real number, e.g., the sampling time. A may be singular. With the exponential of a matrix it is, e.g., possible to compute the solution of a linear system of differential equations

    der(x) = A*x   ->   x(t0 + T) = e^(AT)*x(t0)

The function is called as

     Phi = Matrices.exp(A,T);

       M = Matrices.exp(A);

Algorithmic details:

The algorithm is taken from

H. D. Joos, G. Gruebel:

RASP'91 Regulator Analysis and Synthesis Programs
DLR - Control Systems Group 1991

The following steps are performed to calculate the exponential of A:

1. Matrix A is balanced
   (= is transformed with a diagonal matrix D, such that inv(D)*A*D has a smaller condition as A).
2. The scalar T is divided by a multiple of 2 such that norm( inv(D)*A*D*T/2^k ) < 0.5. Note, that (1) and (2) are implemented such that no round-off errors are introduced.
3. The matrix from (2) is approximated by explicitly performing the Taylor series expansion with a variable number of terms. Truncation occurs if a new term does no longer contribute to the value of Φ from the previous iteration.
4. The resulting matrix is transformed back, by reverting the steps of (2) and (1).

In several sources it is not recommended to use Taylor series expansion to calculate the exponential of a matrix, such as in 'C.B. Moler and C.F. Van Loan: Nineteen dubious ways to compute the exponential of a matrix. SIAM Review 20, pp. 801-836, 1979' or in the documentation of m-file expm2 in Matlab version 6 (http://www.MathWorks.com) where it is stated that 'As a practical numerical method, this is often slow and inaccurate'. These statements are valid for a direct implementation of the Taylor series expansion, but not for the implementation variant used in this function.

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

Inputs

Type	Name	Default	Description
Real	A[:, size(A, 1)]
Real	T	1

Outputs

Type	Name	Description
Real	phi[size(A, 1), size(A, 1)]	= exp(A*T)

Modelica definition

function exp 
  "Compute the exponential of a matrix by adaptive Taylor series expansion with scaling and balancing"

  extends Modelica.Icons.Function;
  input Real A[:, size(A, 1)];
  input Real T=1;
  output Real phi[size(A, 1), size(A, 1)] "= exp(A*T)";

protected 
  parameter Integer nmax=21;
  /*max number of iterations*/
  parameter Integer na=size(A, 1);
  Integer j=1;
  Integer k=0;
  Boolean done=false;
  Real Anorm;
  Real Tscaled=1;
  Real Atransf[na, na];
  Real D[na, na];
  /*D: dummy variable for psi*/
  Real M[na, na];
  /*M: dummy matrix*/
  Real Diag[na];
  /*diagonal transformation matrix for balancing*/

encapsulated function columnNorm 
    "Returns the column norm of a matrix"
  input Real A[:, :] "Input matrix";
  output Real result=0.0 "1-norm of matrix A";
algorithm 
   for i in 1:size(A, 2) loop
      result := max(result, sum(abs(A[:, i])));
   end for;
end columnNorm;

algorithm 
  // balancing of A
  (Diag,Atransf) := balance(A);

  // scaling of T until norm(A)*/(2^k) < 1
  Tscaled := T;
  /*Anorm: column-norm of matrix A*/
  Anorm := columnNorm(Atransf);
  Anorm := Anorm*T;
  while Anorm >= 0.5 loop
    Anorm := Anorm/2;
    Tscaled := Tscaled/2;
    k := k + 1;
  end while;

  // Computation of psi by Taylor-series approximation
  M := identity(na);
  D := M;
  while j < nmax and not done loop
    M := Atransf*M*Tscaled/j;
    //stop if the new element of the series is small
    if columnNorm((D + M) - D) == 0 then
      done := true;
    else
      D := M + D;
      j := j + 1;
    end if;
  end while;

  // re-scaling
  for i in 1:k loop
    D := D*D;
  end for;

  // re-balancing: psi := diagonal(Diag)*D*inv(diagonal(Diag));
  for j in 1:na loop
    for k in 1:na loop
      phi[j, k] := D[j, k]*Diag[j]/Diag[k];
    end for;
  end for;
end exp;

Modelica.Math.Matrices.integralExp

Information

The function uses a Taylor series expansion with Balancing and scaling/squaring to approximate the integral Ψ of the matrix exponential Φ=e^(AT):

                                 AT^2   A^2 * T^3          A^k * T^(k+1)
        Ψ = int(e^(As))ds = IT + ---- + --------- + ... + --------------
                                  2!        3!                (k+1)!

Φ is calculated through Φ = I + A*Ψ, so A may be singular. Γ is simple Ψ*B.

The algorithm runs in the following steps:

1. Balancing
2. Scaling
3. Taylor series expansion
4. Re-scaling
5. Re-Balancing

Balancing put the bad condition of a square matrix A into a diagonal transformation matrix D. This reduce the effort of following calculations. Afterwards the result have to be re-balanced by transformation D*Atransf *inv(D).
Scaling halfen T  k-times, until the norm of A*T is less than 0.5. This garantees minumum rounding errors in the following series expansion. The re-scaling based on the equation  exp(A*2T) = exp(AT)^2. The needed re-scaling formula for psi thus becomes:

         Φ = Φ'*Φ'
   I + A*Ψ = I + 2A*Ψ' + A^2*Ψ'^2
         Ψ = A*Ψ'^2 + 2*Ψ'

where psi' is the scaled result from the series expansion while psi is the re-scaled matrix.

The function is normally used to discretize a state-space system as the zero-order-hold equivalent:

      x(k+1) = Φ*x(k) + Γ*u(k)
        y(k) = C*x(k) + D*u(k)

The zero-order-hold sampling, also known as step-invariant method, gives exact values of the state variables, under the assumption that the control signal u is constant between the sampling instants. Zero-order-hold sampling is discribed in

K. J. Astroem, B. Wittenmark:

Computer Controlled Systems - Theory and Design
Third Edition, p. 32

Syntax:
      (phi,gamma) = Matrices.expIntegral(A,B,T)
                       A,phi: [n,n] square matrices
                     B,gamma: [n,m] input matrix
                           T: scalar, e.g. sampling time

The Algorithm to calculate psi is taken from

H. D. Joos, G. Gruebel:

RASP'91 Regulator Analysis and Synthesis Programs
DLR - Control Systems Group 1991

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

Inputs

Type	Name	Default	Description
Real	A[:, size(A, 1)]
Real	B[size(A, 1), :]
Real	T	1

Outputs

Type	Name	Description
Real	phi[size(A, 1), size(A, 1)]	= exp(A*T)
Real	gamma[size(A, 1), size(B, 2)]	= integral(phi)*B

Modelica definition

function integralExp 
  "Computation of the transition-matrix phi and its integral gamma"

  extends Modelica.Icons.Function;
  input Real A[:, size(A, 1)];
  input Real B[size(A, 1), :];
  input Real T=1;
  output Real phi[size(A, 1), size(A, 1)] "= exp(A*T)";
  output Real gamma[size(A, 1), size(B, 2)] "= integral(phi)*B";
protected 
  parameter Integer nmax=21;
  /*max number of iterations*/
  parameter Integer na=size(A, 1);
  Integer j=2;
  Integer k=0;
  Boolean done=false;
  Real Anorm;
  Real Tscaled=1;
  Real Atransf[na, na];
  Real Psi[na, na];
  /*Psi: dummy variable for psi*/
  Real M[na, na];
  /*M: dummy matrix*/
  Real Diag[na];
  /*diagonal transformation matrix for balancing*/

encapsulated function columnNorm 
    "Returns the column norm of a matrix"
  input Real A[:, :] "Input matrix";
  output Real result=0.0 "1-norm of matrix A";
algorithm 
   for i in 1:size(A, 2) loop
      result := max(result, sum(abs(A[:, i])));
   end for;
end columnNorm;
algorithm 
  // balancing of A
  (Diag,Atransf) := balance(A);

  // scaling of T until norm(A)*/(2^k) < 0.5
  Tscaled := T;
  /*Anorm: column-norm of matrix A*/
  // Anorm := norm(Atransf, 1);
  Anorm := columnNorm(Atransf);
  Anorm := Anorm*T;
  while Anorm >= 0.5 loop
    Anorm := Anorm/2;
    Tscaled := Tscaled/2;
    k := k + 1;
  end while;

  // Computation of psi by Taylor-series approximation
  M := identity(na)*Tscaled;
  Psi := M;
  while j < nmax and not done loop
    M := Atransf*M*Tscaled/j;
    //stop if the new element of the series is small
    // if norm((Psi + M) - Psi, 1) == 0 then
    if columnNorm((Psi + M) - Psi) == 0 then
      done := true;
    else
      Psi := M + Psi;
      j := j + 1;
    end if;
  end while;

  // re-scaling
  for j in 1:k loop
    Psi := Atransf*Psi*Psi + 2*Psi;
  end for;

  // re-balancing: psi := diagonal(Diag)*D*inv(diagonal(Diag));
  for j in 1:na loop
    for k in 1:na loop
      Psi[j, k] := Psi[j, k]*Diag[j]/Diag[k];
    end for;
  end for;
  gamma := Psi*B;
  phi := A*Psi + identity(na);

end integralExp;

Modelica.Math.Matrices.integralExpT

Information

The function calculates the matrices phi,gamma,gamma1 through the equation:

                                 [ A B 0 ]
[phi gamma gamma1] = [I 0 0]*exp([ 0 0 I ]*T)
                                 [ 0 0 0 ]

Syntax:

      (phi,gamma,gamma1) = Matrices.ExpIntegral2(A,B,T)
                     A,phi: [n,n] square matrices
            B,gamma,gamma1: [n,m] matrices
                         T: scalar, e.g. sampling time

The matrices define the discretized first-order-hold equivalent of a state-space system:

      x(k+1) = phi*x(k) + gamma*u(k) + gamma1/T*(u(k+1) - u(k))

K. J. Astroem, B. Wittenmark:

Computer Controlled Systems - Theory and Design
Third Edition, p. 256

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

Inputs

Type	Name	Default	Description
Real	A[:, size(A, 1)]
Real	B[size(A, 1), :]
Real	T	1

Outputs

Type	Name	Description
Real	phi[size(A, 1), size(A, 1)]	= exp(A*T)
Real	gamma[size(A, 1), size(B, 2)]	= integral(phi)*B
Real	gamma1[size(A, 1), size(B, 2)]	= integral((T-t)*exp(A*t))*B

Modelica definition

function integralExpT 
  "Computation of the transition-matrix phi and the integral gamma and gamma1"

  extends Modelica.Icons.Function;
  input Real A[:, size(A, 1)];
  input Real B[size(A, 1), :];
  input Real T=1;
  output Real phi[size(A, 1), size(A, 1)] "= exp(A*T)";
  output Real gamma[size(A, 1), size(B, 2)] "= integral(phi)*B";
  output Real gamma1[size(A, 1), size(B, 2)] "= integral((T-t)*exp(A*t))*B";
protected 
  Integer nmax=200;
  /*max number of iterations*/
  parameter Integer na=size(A, 1);
  parameter Integer nb=size(B, 2);
  Integer j=1;
  Boolean done=false;
  Real F[na + 2*nb, na + 2*nb];

algorithm 
  F := [A, B, zeros(na, nb); zeros(2*nb, na), zeros(2*nb, nb), [identity(nb);
     zeros(nb, nb)]];
  F := exp(F, T);
  phi := F[1:na, 1:na];
  gamma := F[1:na, na + 1:na + nb];
  gamma1 := F[1:na, na + nb + 1:na + 2*nb];

end integralExpT;