## Modelica_LinearSystems2.WorkInProgress

Package for developments and internal models

### Package Content

Name | Description
--- | ---
StateSpace2 | Continuous state space description of a linear, time invariant differential equation system (data + operations)
DiscreteStateSpace | Discrete state space description of a linear, time invariant difference equation system (data + operations)
DiscreteTransferFunction | Discrete transfer function description of a single input, single output system (data + operations)
StateSpace | State space functions
TransferFunction |
Math |
Icons | Icons used in the Linear Systems library
Controller |
Tests | Functions to test the LinearSystems functions

## Modelica_LinearSystems2.WorkInProgress.StateSpace2

Continuous state space description of a linear, time invariant differential equation system (data + operations)

### Information

This record defines a linear time invariant differential equation system in state space form:

        der(x) = A * x + B * u
            y  = C * x + D * u

with

- u - the input vector
- y - the output vector
- x - the state vector
- A,B,C,D - matrices of appropriate dimensions

Extends from Modelica.Icons.Record (Icon for a record).

### Parameters

Type | Name | Default | Description
--- | --- | --- | ---
der(x) = A*x + B*u; y = C*x + D*u
Real | A[:, size(A, 1)] |  |
Real | B[size(A, 1), :] |  |
Real | C[:, size(A, 1)] |  |
Real | D[size(C, 1), size(B, 2)] |  |

### Modelica definition

    record StateSpace2
      "Continuous state space description of a linear, time invariant differential equation system (data + operations)"

      extends Modelica.Icons.Record;

      Real A[:,size(A, 1)];
      Real B[size(A, 1),:];
      Real C[:,size(A, 1)];
      Real D[size(C, 1),size(B, 2)];

    //   String uNames[size(B, 2)]=fill("", size(B, 2))  annotation(Dialog(group="Signal names"));
    //   String yNames[size(C, 1)]=fill("", size(C, 1)) annotation(Dialog(group="Signal names"));
    //   String xNames[size(A, 1)]=fill("", size(A, 1))  annotation(Dialog(group="Signal names"));

    encapsulated operator 'constructor'
        "Default constructors for a StateSpace record"
        import Modelica_LinearSystems2;

      function fromABCDMatrices
          "Default constructor for a StateSpace record"
          import Modelica;
          import Modelica_LinearSystems2;
          import Modelica_LinearSystems2.WorkInProgress.StateSpace2;

        input Real A[:,size(A, 1)];
        input Real B[size(A, 1),:];
        input Real C[:,size(A, 1)];
        input Real D[size(C, 1),size(B, 2)];

          //     input String uNames[size(B, 2)]=fill("", size(B, 2));
          //     input String yNames[size(C, 1)]=fill("", size(C, 1));
          //     input String xNames[size(A, 2)]=fill("", size(A, 2));

         output StateSpace2 result(
          redeclare Real A[size(A, 1),size(A, 2)],
          redeclare Real B[size(B, 1),size(B, 2)],
          redeclare Real C[size(C, 1),size(C, 2)],
          redeclare Real D[size(D, 1),size(D, 2)]);

          //       redeclare String uNames[size(B, 2)],
          //       redeclare String yNames[size(C, 1)],
          //       redeclare String xNames[size(A, 2)]);

      algorithm
        result.A := A;
        result.B := B;
        result.C := C;
        result.D := D;
          //      result.uNames := uNames;
          //      result.yNames := yNames;
          //      result.xNames := xNames;

      end fromABCDMatrices;

      function fromReal
          "Generate a StateSpace data record from a Real value"

          import Modelica;
          import Modelica_LinearSystems2.WorkInProgress.StateSpace2;

        input Real r "Value of Real variable";
        output StateSpace2 ss(
          redeclare Real A[0,0],
          redeclare Real B[0,1],
          redeclare Real C[1,0],
          redeclare Real D[1,1]) "= r";

      algorithm
        ss.D[1, 1] := r;
      end fromReal;

    // algorithm
    //            // this is the contructor algorithm
    //   result.A := A;
    //   result.B := B;
    //   result.C := C;
    //   result.D := D;
    //   //end constructor;

      function fromTransferFunction =
          Modelica_LinearSystems2.TransferFunction.Conversion.toStateSpace;
      function fromZerosAndPoles =
          Modelica_LinearSystems2.ZerosAndPoles.Conversion.toStateSpace;

    end 'constructor';

    encapsulated operator '-'
        "Contains operators for subtraction of state space systems"

      function subtract
          "Subtraction of two state space systems connected in parallel (= inputs are the same, outputs of the two systems are subtracted)"

          import Modelica;
          import Modelica_LinearSystems2.WorkInProgress.StateSpace2;

        input StateSpace2 ss1 "State space system 1";
        input StateSpace2 ss2 "State Space system 2 is subtracted from system 1";
        output StateSpace2 result(
          redeclare Real A[size(ss1.A, 1) + size(ss2.A, 1),size(ss1.A, 2) + size(
            ss2.A, 2)],
          redeclare Real B[size(ss1.B, 1) + size(ss2.B, 1),size(ss1.B, 2)],
          redeclare Real C[size(ss1.C, 1),size(ss1.C, 2) + size(ss2.C, 2)],
          redeclare Real D[size(ss1.D, 1),size(ss1.D, 2)]) "= ss1 - ss2";
        protected
        Integer nx1=size(ss1.A, 1);
        Integer nx2=size(ss2.A, 1);
      algorithm
        result.A := [ss1.A,zeros(nx1, nx2); zeros(nx2, nx1),ss2.A];
        result.B := [ss1.B; ss2.B];
        result.C := [ss1.C,-ss2.C];
        result.D := ss1.D - ss2.D;
      end subtract;

      function negate
          "Unary minus (state space system where the output is multiplied by a gain of -1)"
          import Modelica;
          import StateSpace = Modelica_LinearSystems2.WorkInProgress.StateSpace2;

        input StateSpace ss;
        output StateSpace result(
          redeclare Real A[size(ss.A, 1),size(ss.A, 2)],
          redeclare Real B[size(ss.B, 1),size(ss.B, 2)],
          redeclare Real C[size(ss.C, 1),size(ss.C, 2)],
          redeclare Real D[size(ss.D, 1),size(ss.D, 2)]) "= -ss";
      algorithm
        result.A := ss.A;
        result.B := ss.B;
        result.C := -ss.C;
        result.D := -ss.D;
      end negate;
    end '-';

    encapsulated operator function '+'
        "Parallel connection of two state space systems (= inputs are the same, outputs of the two systems are added)"
        import StateSpace = Modelica_LinearSystems2.WorkInProgress.StateSpace2;

        input StateSpace ss1 "System 1";
        input StateSpace ss2 "System 2 is added in parallel to system 1";
        output StateSpace result(
          redeclare Real A[size(ss1.A, 1) + size(ss2.A, 1),size(ss1.A, 2) + size(
            ss2.A, 2)],
          redeclare Real B[size(ss1.B, 1) + size(ss2.B, 1),size(ss1.B, 2)],
          redeclare Real C[size(ss1.C, 1),size(ss1.C, 2) + size(ss2.C, 2)],
          redeclare Real D[size(ss1.D, 1),size(ss1.D, 2)]) "= ss1 + ss2";
      protected
        Integer nx1=size(ss1.A, 1);
        Integer nx2=size(ss2.A, 1);
    algorithm
        result.A := [ss1.A,zeros(nx1, nx2); zeros(nx2, nx1),ss2.A];
        result.B := [ss1.B; ss2.B];
        result.C := [ss1.C,ss2.C];
        result.D := ss1.D + ss2.D;
      //end '+';
    //  end add;

    end '+';

    encapsulated operator function '*'
        "Series connection of two state space systems"
        import StateSpace = Modelica_LinearSystems2.WorkInProgress.StateSpace2;

        input StateSpace ss1 "System 1";
        input StateSpace ss2 "System 2";
        output StateSpace result(
          redeclare Real A[size(ss1.A, 1) + size(ss2.A, 1),size(ss1.A, 2) + size(
            ss2.A, 2)],
          redeclare Real B[size(ss1.B, 1) + size(ss2.B, 1),size(ss2.B, 2)],
          redeclare Real C[size(ss1.C, 1),size(ss1.C, 2) + size(ss2.C, 2)],
          redeclare Real D[size(ss1.D, 1),size(ss2.D, 2)])
          "y = G(s)*u = G(ss1)*G(ss2)*u";
      protected
        Integer nx1=size(ss1.A, 1);
        Integer nx2=size(ss2.A, 1);
    algorithm
        result.A := [ss1.A,ss1.B*ss2.C; zeros(nx2, nx1),ss2.A];
        result.B := [ss1.B*ss2.D; ss2.B];
        result.C := [ss1.C,ss1.D*ss2.C];
        result.D := ss1.D*ss2.D;
    //  end multiply;

    end '*';

    encapsulated operator function '=='
        "Check whether two linear systems have identical matrices"
        import Modelica.Math.Matrices.isEqual;
        import StateSpace = Modelica_LinearSystems2.WorkInProgress.StateSpace2;

        input StateSpace ss1 "System 1";
        input StateSpace ss2 "System 2";
        input Real eps(min=0) = 0
          "Two elements e1 and e2 of the two systems are identical if abs(e1-e2) <= eps";
        output Boolean same "=true, if the two systems are identical";
    algorithm
        same := isEqual(
              ss1.A,
              ss2.A,
              eps) and isEqual(
              ss1.B,
              ss2.B,
              eps) and isEqual(
              ss1.C,
              ss2.C,
              eps) and isEqual(
              ss1.D,
              ss2.D,
              eps);
    end '==';

    encapsulated package Import
        "Utilitiy functions to import StaeSpace representations"

      encapsulated function fromFile
          "Read a StateSpace data record from mat-file"

          import Modelica;
          import Modelica_LinearSystems2.StateSpace;
          import Modelica_LinearSystems2.WorkInProgress.StateSpace2;
          import Modelica_LinearSystems2;

        input String fileName="dslin.mat"
            "Name of the state space system data file";
        input String matrixName="ABCD" "Name of the state space system matrix";
        protected
        input Integer xuy[3]=Modelica_LinearSystems2.StateSpace.Internal.readSystemDimension(fileName, matrixName);
        input Integer nx=xuy[1];
        input Integer nu=xuy[2];
        input Integer ny=xuy[3];

        public
        output StateSpace2 result(
          redeclare Real A[nx,nx],
          redeclare Real B[nx,nu],
          redeclare Real C[ny,nx],
          redeclare Real D[ny,nu]) "= model linearized at initial point";

        protected
        Real ABCD[nx + ny,nx + nu]=Modelica_LinearSystems2.Internal.Streams.readMatrixInternal(
              fileName,
              matrixName,
              nx + ny,
              nx + nu);

      algorithm
        result.A := ABCD[1:nx, 1:nx];
        result.B := ABCD[1:nx, nx + 1:nx + nu];
        result.C := ABCD[nx + 1:nx + ny, 1:nx];
        result.D := ABCD[nx + 1:nx + ny, nx + 1:nx + nu];
        Modelica.Utilities.Streams.print("StateSpace record loaded from file: \""
           + fileName + "\"");

      end fromFile;

    end Import;

    end StateSpace2;

## Modelica_LinearSystems2.WorkInProgress.DiscreteStateSpace

Discrete state space description of a linear, time invariant difference equation system (data + operations)

### Information

This record defines a linear time invariant difference equation system in state space form:

         x(Ts*(k+1)) = A * x(Ts*k) + B * u(Ts*k)
         y(Ts*k)     = C * x(Ts*k) + D * u(Ts*k)
         x_continuous(Ts*k) = x(Ts*k) + B2 * u(Ts*k)

with

- Ts - the sample time
- k - the index of the actual sample instance (k=0,1,2,3,...)
- t - the time
- u(t) - the input vector,
- y(t) - the output vector,
- x(t) - the discrete state vector (x(t=Ts*0) is the initial state),
- x_continuous(t) - the state vector of the continuous system from which the discrete block has been derived (details see below),
- A,B,C,D,B2 - matrices of appropriate dimensions.

A discrete system is usually derived by discretization from a continuous block, e.g., by function LinearSystems.DiscreteStateSpace.fromStateSpace. If the discretization method, e.g., the trapezoidal method, accesses actual and past values of the input u (e.g. u(Ts*k), u(Ts*(k-1), u(Ts*(k-2))), a state transformation is needed to get the difference equation above where only the actual value u(Ts*k) is accessed.

If the original continuous state vector should be computed from the sampled data system above, the matrices of this transformation have to be known. For simplicity and efficiency, here only the specific transformation used by function LinearSystems.DiscreteStateSpace.fromStateSpace is stored in the data record of the discrete system via matrix B2. Therefore, the state vector of the underlying continuous system can be calculated by adding the term B2*u to the state vector of the discretized system.

In Modelica notation, the difference equation above can be implemented as:

         when {initial(), sample(Ts,Ts)} then
            new_x = A * x + B * u;
                y = C * x + D * u;
                x = pre(new_x);
            x_continuous = x + B2 * u;
         end when;

Since no "next value" operator is available in Modelica, an auxiliary variable new_x stores the value of x for the next sampling instant. The relationship between new_x and x is defined via equation "x = pre(new_x)".

The body of the when-clause is active during initialization and at the next sample instant t=Ts. Note, the when-equation is not active after the initialization at t=0 (due to sample(Ts,Ts)), since the state x of the initialization has to be used also at t=0.

In library Blocks.Controller additional equations are added for the initialization to uniquely compute the initial vector x:

      initial equation
         if init == InitialState then
            x = x_start;
         elseif init == SteadyState then
            x = new_x;
         end if;

Optionally, x is set to a given start vector x_start. As default initialization, the equation "x = new_x" is added that defines steady state initialization for the discrete system. As a consequence, the output y(Ts*k), k=0,1,2,.., remains constant after this initialization, provided the input vector u(Ts*k) remains constant.

Extends from Modelica.Icons.Record (Icon for a record).

### Parameters

Type | Name | Default | Description
--- | --- | --- | ---
new_x = A*x + B*u; y = C*x + D*u; x_cont = x + B2*u
Real | A[:, size(A, 1)] |  |
Real | B[size(A, 1), :] |  |
Real | C[:, size(A, 1)] |  |
Real | D[size(C, 1), size(B, 2)] |  |
Data used to construct discrete from continuous system
Time | Ts | 1 | Sample time [s]
Real | B2[size(B, 1), size(B, 2)] | fill(0, size(B, 1), size(B, ... | Reconstruct continuous state
Method | method | Modelica_LinearSystems2.Type... | Discretization method

### Modelica definition

    record DiscreteStateSpace
      "Discrete state space description of a linear, time invariant difference equation system (data + operations)"
      extends Modelica.Icons.Record;

      Real A[:,size(A, 1)];
      Real B[size(A, 1),:];
      Real C[:,size(A, 1)];
      Real D[size(C, 1),size(B, 2)];

      Modelica.SIunits.Time Ts=1 "Sample time";
      Real B2[size(B, 1),size(B, 2)]=fill(0,size(B,1),size(B,2))
        "Reconstruct continuous state";
      Modelica_LinearSystems2.Types.Method method=Modelica_LinearSystems2.Types.Method.Trapezoidal
        "Discretization method";

      encapsulated operator 'constructor'
        function fromMatrices
          "Default constructor for a DiscreteStateSpace record"
          import Modelica;
          import Modelica_LinearSystems2;
          import Modelica_LinearSystems2.WorkInProgress.DiscreteStateSpace;

          input Real A[:,size(A, 1)];
          input Real B[size(A, 1),:];
          input Real C[:,size(A, 1)];
          input Real D[size(C, 1),size(B, 2)];

          input Modelica.SIunits.Time Ts=1 "Sample time";
          input Real B2[:,:]=zeros(size(B, 1), size(B, 2));
          input Modelica_LinearSystems2.Types.Method method=Modelica_LinearSystems2.Types.Method.Trapezoidal
            "Discretization methodDiscretization method";
          output DiscreteStateSpace result(
            redeclare Real A[size(A, 1),size(A, 2)],
            redeclare Real B[size(B, 1),size(B, 2)],
            redeclare Real C[size(C, 1),size(C, 2)],
            redeclare Real D[size(D, 1),size(D, 2)],
            redeclare Real B2[size(B2, 1),size(B2, 2)]);
        algorithm
          result.A := A;
          result.B := B;
          result.C := C;
          result.D := D;

          result.B2 := B2;
          result.Ts := Ts;
          result.method := method;
        end fromMatrices;

        function fromStateSpace
          "Transform a continuous into a discrete linear state space system"
          import Modelica;
          import Modelica_LinearSystems2;
          import Modelica_LinearSystems2.Types.Method;
          import Modelica_LinearSystems2.Math.Matrices.LU_solve2;

          input Modelica_LinearSystems2.StateSpace sc
            "Continuous linear state space system";
          input Modelica.SIunits.Time Ts "Sample time";
          input Modelica_LinearSystems2.Types.Method method=Modelica_LinearSystems2.Types.Method.Trapezoidal
            "Discretization method";
          output Modelica_LinearSystems2.WorkInProgress.DiscreteStateSpace sd(
            redeclare Real A[size(sc.A, 1),size(sc.A, 2)],
            redeclare Real B[size(sc.B, 1),size(sc.B, 2)],
            redeclare Real C[size(sc.C, 1),size(sc.C, 2)],
            redeclare Real D[size(sc.D, 1),size(sc.D, 2)],
            redeclare Real B2[size(sc.B, 1),size(sc.B, 2)])
            "Discrete state space system";

        protected
          Integer nx=size(sc.A, 1) "Number of states";
          Integer nu=size(sc.B, 2) "Number of input signals";
          Real LU[nx,nx] "LU decomposition";
          Integer pivots[nx] "Pivots of LU decomposition";
        algorithm
          sd.Ts := Ts;
          sd.method := method;

          if method == Method.ExplicitEuler then
                /*  der_x = A*x + B*u
                 x = pre(x) + Ts*pre(der_x)
         */
            sd.A := identity(nx) + Ts*sc.A;
            sd.B := Ts*sc.B;
            sd.C := sc.C;
            sd.D := sc.D;
            sd.B2 := zeros(nx, nu);

          elseif method == Method.ImplicitEuler then
                /*  der_x = A*x + B*u
                 x = pre(x) + Ts*der_x
         */
            (LU,pivots) := Modelica_LinearSystems2.Math.Matrices.LU(identity(nx) -
              Ts*sc.A);
            sd.B2 := LU_solve2(
                  LU,
                  pivots,
                  Ts*sc.B);
            sd.A := LU_solve2(
                  LU,
                  pivots,
                  identity(nx));
            sd.B := sd.A*sd.B2;
            sd.C := sc.C;
            sd.D := sd.C*sd.B2 + sc.D;

          elseif method == Method.Trapezoidal then
                /*  der_x = A*x + B*u
                 x = pre_x + (Ts/2)*(pre_der_x + der_x);
         */
            (LU,pivots) := Modelica_LinearSystems2.Math.Matrices.LU(identity(nx) -
              (Ts/2)*sc.A);
            sd.B2 := LU_solve2(
                  LU,
                  pivots,
                  (Ts/2)*sc.B);
            sd.A := LU_solve2(
                  LU,
                  pivots,
                  identity(nx) + (Ts/2)*sc.A);
            sd.B := sd.A*sd.B2 + sd.B2;
            sd.C := sc.C;
            sd.D := sd.C*sd.B2 + sc.D;

          elseif method == Method.StepExact then
               /* x = phi*pre(x) + gamma*pre(u);
           y = C*x + D*u
        */
            (sd.A,sd.B) := Modelica.Math.Matrices.integralExp(
                  sc.A,
                  sc.B,
                  Ts);
            sd.C := sc.C;
            sd.D := sc.D;
            sd.B2 := zeros(nx, nu);

          elseif method == Method.RampExact then
               /* x = phi*pre(x) + gamma*pre(u) + gamma1/Ts*(u - pre_u);
            -> x = phi*pre(x) + (gamma - gamma1/Ts)*pre(u) + gamma1/Ts*u;
      z = x - gamma1/Ts*u
          leads to
            z = phi*pre(z) + phi*gamma1/Ts*pre(u) +gamma*pre(u) - gamma1/Ts*pre(u)
           y = C*x + D*u
          -> y = C*z + (D + C*gamma1/Ts)*u
          x = z + gamma1/Ts*u

        */
            (sd.A,sd.B,sd.B2) := Modelica.Math.Matrices.integralExpT(
                  sc.A,
                  sc.B,
                  Ts);
            sd.B2 := sd.B2/Ts;
            sd.B := sd.A*sd.B2 + sd.B - sd.B2;
            sd.C := sc.C;
            sd.D := sd.C*sd.B2 + sc.D;

          elseif method == Method.ImpulseExact then
               /* x = phi*pre(x) + phi*B*u;
            y = C*x
            (u = [1,0,1,1,0..,0])
          Limitations: The infinit impulses at t = kT is ignored in the mapping
        */

            sd.A := Modelica.Math.Matrices.exp(sc.A, Ts);
            sd.B := sd.A*sc.B;
            sd.C := sc.C;
            sd.D := sc.C*sc.B;
            sd.B2 := sc.B;

          else
            assert(false, "Argument method (= " + String(method) +
              ") of makeDiscrete is wrong.");
          end if;
        end fromStateSpace;

        encapsulated function fromMatrices2
          "Transform a continuous into a discrete linear state space system"
          import Modelica;
          import Modelica_LinearSystems2;
          import Modelica_LinearSystems2.Types.Method;
          import Modelica_LinearSystems2.Math.Matrices.LU_solve2;

          input Real A[:,size(A, 1)];
          input Real B[size(A, 1),:];
          input Real C[:,size(A, 1)];
          input Real D[size(C, 1),size(B, 2)];
          input Modelica.SIunits.Time Ts "Sample time";
          input Modelica_LinearSystems2.Types.Method method=Modelica_LinearSystems2.Types.Method.Trapezoidal
            "Discretization method";
        //  input Modelica_LinearSystems2.Types method=Modelica_LinearSystems2.Types.Method.Trapezoidal
          output Modelica_LinearSystems2.WorkInProgress.DiscreteStateSpace sd(
            redeclare Real A[size(A, 1),size(A, 2)],
            redeclare Real B[size(B, 1),size(B, 2)],
            redeclare Real C[size(C, 1),size(C, 2)],
            redeclare Real D[size(D, 1),size(D, 2)],
            redeclare Real B2[size(B, 1),size(B, 2)]) "Discrete state space system";

        protected
          Integer nx=size(A, 1) "Number of states";
          Integer nu=size(B, 2) "Number of input signals";
          Real LU[nx,nx] "LU decomposition";
          Integer pivots[nx] "Pivots of LU decomposition";
        algorithm
          sd.Ts := Ts;
          sd.method := method;

          if method == Method.ExplicitEuler then
                /*  der_x = A*x + B*u
                 x = pre(x) + Ts*pre(der_x)
         */
            sd.A := identity(nx) + Ts*A;
            sd.B := Ts*B;
            sd.C := C;
            sd.D := D;
            sd.B2 := zeros(nx, nu);

          elseif method == Method.ImplicitEuler then
                /*  der_x = A*x + B*u
                 x = pre(x) + Ts*der_x
         */
            (LU,pivots) := Modelica_LinearSystems2.Math.Matrices.LU(identity(nx) -
              Ts*A);
            sd.B2 := LU_solve2(
                  LU,
                  pivots,
                  Ts*B);
            sd.A := LU_solve2(
                  LU,
                  pivots,
                  identity(nx));
            sd.B := sd.A*sd.B2;
            sd.C := C;
            sd.D := sd.C*sd.B2 + D;

          elseif method == Method.Trapezoidal then
                /*  der_x = A*x + B*u
                 x = pre_x + (Ts/2)*(pre_der_x + der_x);
         */
            (LU,pivots) := Modelica_LinearSystems2.Math.Matrices.LU(identity(nx) -
              (Ts/2)*A);
            sd.B2 := LU_solve2(
                  LU,
                  pivots,
                  (Ts/2)*B);
            sd.A := LU_solve2(
                  LU,
                  pivots,
                  identity(nx) + (Ts/2)*A);
            sd.B := sd.A*sd.B2 + sd.B2;
            sd.C := C;
            sd.D := sd.C*sd.B2 + D;

          elseif method == Method.StepExact then
               /* x = phi*pre(x) + gamma*pre(u);
           y = C*x + D*u
        */
            (sd.A,sd.B) := Modelica.Math.Matrices.integralExp(
                  A,
                  B,
                  Ts);
            sd.C := C;
            sd.D := D;
            sd.B2 := zeros(nx, nu);

          elseif method == Method.RampExact then
               /* x = phi*pre(x) + gamma*pre(u) + gamma1/Ts*(u - pre_u);
            -> x = phi*pre(x) + (gamma - gamma1/Ts)*pre(u) + gamma1/Ts*u;
      z = x - gamma1/Ts*u
          leads to
            z = phi*pre(z) + phi*gamma1/Ts*pre(u) +gamma*pre(u) - gamma1/Ts*pre(u)
           y = C*x + D*u
          -> y = C*z + (D + C*gamma1/Ts)*u
          x = z + gamma1/Ts*u

        */
            (sd.A,sd.B,sd.B2) := Modelica.Math.Matrices.integralExpT(
                  A,
                  B,
                  Ts);
            sd.B2 := sd.B2/Ts;
            sd.B := sd.A*sd.B2 + sd.B - sd.B2;
            sd.C := C;
            sd.D := sd.C*sd.B2 + D;

          elseif method == Method.ImpulseExact then
               /* x = phi*pre(x) + phi*B*u;
            y = C*x
            (u = [1,0,1,1,0..,0])
          Limitations: The infinit impulses at t = kT is ignored in the mapping
        */

            sd.A := Modelica.Math.Matrices.exp(A, Ts);
            sd.B := sd.A*B;
            sd.C := C;
            sd.D := C*B;
            sd.B2 := B;

          else
            assert(false, "Argument method (= " + String(method) +
              ") of makeDiscrete is wrong.");
          end if;
        end fromMatrices2;
      end 'constructor';

      encapsulated function timeResponse
        "Compute time response of DiscreteStateSpace system"
        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.WorkInProgress.DiscreteStateSpace;

        input DiscreteStateSpace sd "Linear system in discrete state space form";
        input Real u[:,size(sd.B, 2)]=ones(3, size(sd.B, 2))
          "System input (dimension: (input samples) x (number of inputs))";
        input Real x0[size(sd.A, 1)]=zeros(size(sd.A, 1)) "Initial system state";
        output Real y[size(u, 1),size(sd.C, 1)]
          "System response (dimension: (input samples) x (number of outputs))";
        output Real x_continuous[size(u, 1),size(sd.A, 1)]
          "State trajectories (dimension: (input samples) x (number of states)";

      protected
        Integer samples=size(u, 1);
        Integer i;
        Real new_x[size(sd.A, 1)];
        Real x[size(sd.A, 1)]=x0;

      algorithm
        for i in 1:samples loop
          new_x := sd.A*x + sd.B*u[i, :];
          y[i, :] := sd.C*x + sd.D*u[i, :];
          x_continuous[i, :] := x + sd.B2*u[i, :];
          x := new_x;
        end for;
      end timeResponse;

      encapsulated function initialResponse
        "Compute initial response of DiscreteStateSpace system"
        import Modelica;
        import Modelica_LinearSystems2;
        import Modelica_LinearSystems2.WorkInProgress.DiscreteStateSpace;

        input DiscreteStateSpace sd "Linear system in discrete state space form";
        input Real x0[size(sd.A, 1)]=zeros(size(sd.A, 1)) "Initial system state";
        input Integer samples "Number of samples";
        output Real y[samples,size(sd.C, 1)]
          "System response (dimension: (input samples) x (number of outputs))";
        output Real x_continuous[samples,size(sd.A, 1)]
          "State trajectories (dimension: (input samples) x (number of states)";

      protected
        Integer i;
        Real new_x[size(sd.A, 1)];
        Real x[size(sd.A, 1)]=x0;

      algorithm
        for i in 1:samples loop
          new_x := sd.A*x;
          y[i, :] := sd.C*x;
          x_continuous[i, :] := x;
          x := new_x;
        end for;
      end initialResponse;

    encapsulated package Conversion
        "Conversion functions from DiscreteStateSpace into DiscreteTransferFunction"
    function toDiscreteTransferFunction
          "Generate a transfer function from a SISO state space representation"

          import Modelica;
          import Modelica_LinearSystems2;
          import Modelica_LinearSystems2.TransferFunction;
          import Modelica_LinearSystems2.StateSpace;
          import Modelica_LinearSystems2.WorkInProgress.DiscreteTransferFunction;
          import Modelica_LinearSystems2.WorkInProgress.DiscreteStateSpace;

      input DiscreteStateSpace dss "DiscreteStateSpace object";

      output Modelica_LinearSystems2.WorkInProgress.DiscreteTransferFunction dtf
            "DiscreteTransferFunction object";

        protected
      StateSpace ss = StateSpace(dss.A,dss.B,dss.C,dss.D);
      TransferFunction tf = StateSpace.Conversion.toTransferFunction(ss);

    algorithm
      dtf := DiscreteTransferFunction(n=tf.n, d=tf.d, Ts=dss.Ts, method=dss.method);

    end toDiscreteTransferFunction;
    end Conversion;

    encapsulated package Plot
        "Functions to plot discrete state space system responses"
    encapsulated function bodeSISO
          "Plot bode plot of the corresponding discrete transfer function"
          import Modelica;
          import Modelica_LinearSystems2;
          import Modelica_LinearSystems2.WorkInProgress.DiscreteStateSpace;
          import Modelica_LinearSystems2.WorkInProgress.DiscreteTransferFunction;

      input DiscreteStateSpace dss "discrete state space system";
      input Integer iu=1 "index of input";
      input Integer iy=1 "index of output";
      input Integer nPoints(min=2) = 200 "Number of points";
      input Boolean autoRange=true
            "= true, if abszissa range is automatically determined";
      input Modelica.SIunits.Frequency f_min=0.1
            "Minimum frequency value, if autoRange = false";
      input Modelica.SIunits.Frequency f_max=10
            "Maximum frequency value, if autoRange = false";

      input Boolean magnitude=true "= true, to plot the magnitude of dtf";
      input Boolean phase=true "= true, to plot the pase of tf";

      extends Modelica_LinearSystems2.Internal.PartialPlotFunction(defaultDiagram=
            Modelica_LinearSystems2.Internal.DefaultDiagramBodePlot());

        protected
      Modelica_LinearSystems2.WorkInProgress.DiscreteTransferFunction dtf
            "Transfer functions to be plotted";
      DiscreteStateSpace dss_siso(
        redeclare Real A[size(dss.A, 1),size(dss.A, 2)],
        redeclare Real B[size(dss.B, 1),1],
        redeclare Real C[1,size(dss.C, 2)],
        redeclare Real D[1,1],
        redeclare Real B2[size(dss.B, 1),1]);

    algorithm
      assert(iu <= size(dss.B, 2) and iu > 0, "index for input is " + String(iu) + " which is not in [1, "
         + String(size(dss.B, 2)) + "].");
      assert(iy <= size(dss.C, 1) and iy > 0, "index for output is " + String(iy) + " which is not in [1, "
         + String(size(dss.C, 1)) + "].");
      dss_siso := DiscreteStateSpace(
        A=dss.A,
        B=matrix(dss.B[:, iu]),
        C=transpose(matrix(dss.C[iy, :])),
        D=matrix(dss.D[iy, iu]),
        B2=matrix(dss.B2[:, iu]),
        Ts=dss.Ts,
        method=dss.method);
      dtf := DiscreteStateSpace.Conversion.toDiscreteTransferFunction(dss_siso);

      DiscreteTransferFunction.Plot.bode(
        dtf,
        nPoints,
        autoRange,
        f_min,
        f_max,
        magnitude,
        phase,
        defaultDiagram=defaultDiagram,
        device=device);

    end bodeSISO;

    end Plot;
    end DiscreteStateSpace;

## Modelica_LinearSystems2.WorkInProgress.DiscreteTransferFunction

Discrete transfer function description of a single input, single output system (data + operations)

### Information

Extends from Modelica.Icons.Record (Icon for a record).

### Parameters

Type | Name | Default | Description
--- | --- | --- | ---
y = n*{s^m, ... , s, 1} / (d*{s^r, ... , s, 1}) * u
Real | n[:] |  | Coefficients of numerator polynomial (in descending order)
Real | d[:] |  | Coefficients of denominator polynomial (in descending order)
Data used to construct discrete from continuous system
Time | Ts |  | Sample time [s]
Method | method | Modelica_LinearSystems2.Type... | Discretization method
Signal names
String | uName | "u" | Name of input signal
String | yName | "y" | Name of output signal

### Modelica definition

    record DiscreteTransferFunction
      "Discrete transfer function description of a single input, single output system (data + operations)"

      extends Modelica.Icons.Record;
      import Modelica_LinearSystems2.Math.Polynomial;
      import Modelica_LinearSystems2;

      Real n[:] "Coefficients of numerator polynomial (in descending order)";
      Real d[:] "Coefficients of denominator polynomial (in descending order)";

      Modelica.SIunits.Time Ts "Sample time";

      Modelica_LinearSystems2.Types.Method method=Modelica_LinearSystems2.Types.Method.Trapezoidal
        "Discretization method";

      String uName="u" "Name of input signal";
      String yName="y" "Name of output signal";

    /* If the numerator polynomial has no coefficients, the transfer function
       is zero. The denominator polynomial must always have at
       least one coefficient, such as {1}
    */

      encapsulated operator 'constructor'
        "Default constructor for a transfer function"
        import Modelica;
        import Modelica_LinearSystems2.TransferFunction;

        function fromReal
          "Generate a DiscreteTransferFunction data record from a Real value"
          import Modelica;
          import Modelica_LinearSystems2.WorkInProgress.DiscreteTransferFunction;

          input Real r "Value of Real variable";
          input String uName="" "input name";
          input String yName="" "output name";
          output DiscreteTransferFunction dtf(n={r}, d={1});

        algorithm
          dtf.uName := uName;
          dtf.yName := yName;
        end fromReal;

          encapsulated function fromArrays
          "Generate a DiscreteTransferFunction data record from numerator and denominator array"
          import Modelica;
          import Modelica_LinearSystems2.WorkInProgress.DiscreteTransferFunction;
          import Modelica_LinearSystems2;

               input Real n[:] "Coefficients of numerator polynomial";
               input Real d[:] "Coefficients of denominator polynomial";
               input Modelica.SIunits.Time Ts "Sample time";
               input Modelica_LinearSystems2.Types.Method method=Modelica_LinearSystems2.Types.Method.Trapezoidal
            "Discretization method";

               input String uName = "" "input name";
               input String yName = "" "output name";

               output DiscreteTransferFunction dtf(redeclare Real n[size(n, 1)], redeclare
              Real d[                                                                             size(d, 1)])
            "Transfer function";

          algorithm
                     //this is the constructor algorithm
               assert(size(d, 1) > 0, "Input denominator d must have at least one element, however\n"
                  + "d is an empty vector");
               dtf.n := n;
               dtf.d := d;
               dtf.Ts := Ts;
               dtf.method := method;
               dtf.uName := uName;
               dtf.yName := yName;
          end fromArrays;

        function fromPolynomials
          "Generate a DiscreteTransferFunction data record from a numerator and denominator polynomial"
          import Modelica;
          import Modelica_LinearSystems2.WorkInProgress.DiscreteTransferFunction;
          import Modelica_LinearSystems2.Math.Polynomial;

          input Polynomial n "Numerator polynomial";
          input Polynomial d "Denominator polynomial";
          input String uName="" "input name";
          input String yName="" "output name";
          output DiscreteTransferFunction dtf(n=n.c, d=d.c,uName=uName, yName=yName);

        algorithm

        end fromPolynomials;

        function fromTransferFunction
          "Generate a DiscreteTransferFunction data record from a continuous Transfer function"
          import Modelica;
          import Modelica_LinearSystems2;
          import Modelica_LinearSystems2.WorkInProgress.DiscreteTransferFunction;
          import Modelica_LinearSystems2.TransferFunction;
          import Modelica_LinearSystems2.StateSpace;
          import Modelica_LinearSystems2.WorkInProgress.DiscreteStateSpace;

          input TransferFunction tf "continuous transfer function";
          input Modelica.SIunits.Time Ts "Sample time";

          input Modelica_LinearSystems2.Types.Method method=Modelica_LinearSystems2.Types.Method.Trapezoidal
            "Discretization method";

          output DiscreteTransferFunction dtf;
        protected
          StateSpace ss = StateSpace(tf);
          Modelica_LinearSystems2.WorkInProgress.DiscreteStateSpace dss=
                                   Modelica_LinearSystems2.WorkInProgress.DiscreteStateSpace(
                                                      ss,Ts,method);

        algorithm
          dtf := DiscreteStateSpace.Conversion.toDiscreteTransferFunction(dss);
        end fromTransferFunction;
      end 'constructor';

        encapsulated operator function 'String'
        "Transform TransferFunction into a String representation"
        import Modelica_LinearSystems2.Math.Polynomial;
        import Modelica_LinearSystems2.WorkInProgress.DiscreteTransferFunction;

        input DiscreteTransferFunction dtf
          "Discrete transfer function to be transformed in a String representation";
        input Integer significantDigits=6
          "Number of significant digits that are shown";
        input String name="z" "Independent variable name used for printing";
        output String z="";
      protected
        Integer n_num=size(dtf.n, 1) - 1;
        Integer n_den=size(dtf.d, 1) - 1;
        Boolean numParenthesis;
        algorithm
        if n_num == -1 then
          z := "0";
        else
          numParenthesis := n_num > 0 and not (n_den == 0 and dtf.d[1] == 1);
          if numParenthesis then
            z := "(";
          end if;
           z := z + String(
                Polynomial(dtf.n),
                significantDigits,
                name);

          if numParenthesis then
            z := z + ")";
          end if;
        if n_den > 0 or dtf.d[1] <> 1 then
            if n_den > 0 then
              z := z + "/(";
            else
              z := z + "/";
            end if;

            z := z + String(
                  Polynomial(dtf.d),
                  significantDigits,
                  name);

            if n_den > 0 then
              z := z + ")";
            end if;
          end if;
        end if;
        end 'String';

    encapsulated package Plot
        "Functions to plot state space system responses"

    encapsulated function bode "Plot transfer function as bode plot"
          import Modelica;
          import Modelica.Utilities.Strings;
          import Modelica_LinearSystems2;
          import Modelica_LinearSystems2.Internal;
          import Modelica_LinearSystems2.TransferFunction;
          import Modelica_LinearSystems2.WorkInProgress.DiscreteTransferFunction;
          import Modelica_LinearSystems2.Math.Complex;

          import Modelica_LinearSystems2.Utilities.Plot;

          import SI = Modelica.SIunits;

      input DiscreteTransferFunction dtf "DiscreteTransfer function to be plotted";
      input Integer nPoints(min=2) = 200 "Number of points";
      input Boolean autoRange=true
            "= true, if abszissa range is automatically determined";
      input Modelica.SIunits.Frequency f_min(min=0) = 0.1
            "Minimum frequency value, if autoRange = false";
      input Modelica.SIunits.Frequency f_max(min=0) = 10
            "Maximum frequency value, if autoRange = false";

      input Boolean magnitude=true "= true, to plot the magnitude of tf";
      input Boolean phase=true "= true, to plot the pase of tf";

      extends Modelica_LinearSystems2.Internal.PartialPlotFunction(defaultDiagram=
            Modelica_LinearSystems2.Internal.DefaultDiagramBodePlot(heading="Bode plot of  dtf = "
             + String(dtf)));

        protected
      SI.AngularVelocity w[nPoints];
      Complex z[nPoints];
      SI.Frequency f[nPoints];
      Modelica.SIunits.Conversions.NonSIunits.Angle_deg phi[nPoints];
      Real A[nPoints];
      Boolean OK;
      Complex c;
      Modelica.SIunits.Angle phi_old;
      Complex numZeros[:];
      Complex denZeros[:];
      Complex numZerosZ[:];
      Complex denZerosZ[:];
      TransferFunction tf=TransferFunction(n=dtf.n, d=dtf.d);

      Plot.Records.Curve curves[2];
      Integer i;
      Plot.Records.Diagram diagram2[2];

    algorithm
            // Determine frequency vector f
      if autoRange then
        (numZerosZ,denZerosZ) := TransferFunction.Analysis.zerosAndPoles(tf);
      else
        numZerosZ := fill(Complex(0), 0);
        denZerosZ := fill(Complex(0), 0);
      end if;

      numZeros := fill(Complex(0),0);
     // numZeros := fill(Complex(0),size(numZerosZ,1));
      // for i in 1:size(numZerosZ,1) loop
      //   numZeros[i] := Complex.log(numZerosZ[i])/dtf.Ts;
      // end for;

      denZeros := fill(Complex(0),size(denZerosZ,1));
      for i in 1:size(denZerosZ,1) loop
        denZeros[i] := Complex.log(denZerosZ[i])/dtf.Ts;
      end for;

      f := Internal.frequencyVector(
            nPoints,
            autoRange,
            f_min,
            f_max,
            numZeros,
            denZeros);

      // Compute magnitude/phase at the frequency points
      phi_old := 0.0;
      for i in 1:nPoints loop
        w[i] := Modelica.SIunits.Conversions.from_Hz(f[i]);
        z[i] := Complex.exp(Complex(0,w[i]*dtf.Ts));
        c := TransferFunction.Analysis.evaluate(
              tf,
              z[i],
              1e-10);
        A[i] := Complex.'abs'(c);
        phi_old := Complex.arg(c, phi_old);
        phi[i] := Modelica.SIunits.Conversions.to_deg(phi_old);

      end for;

     // Plot computed frequency response
     diagram2 := fill(defaultDiagram, 2);
      i := 0;
      if magnitude then
        i := i + 1;
        curves[i] := Plot.Records.Curve(
              x=f,
              y=A,
              autoLine=true);
        diagram2[i].curve := {curves[i]};
        diagram2[i].yLabel := "magnitude";
        if phase then
           diagram2[i].xLabel:="";
        end if;
      end if;

      if phase then
        i := i + 1;
        curves[i] := Plot.Records.Curve(
              x=f,
              y=phi,
              autoLine=true);
        diagram2[i].curve := {curves[i]};
        diagram2[i].yLabel := "phase [deg]";
        diagram2[i].logY := false;
        if magnitude then
          diagram2[i].heading:="";
       end if;
      end if;

        if magnitude and phase then
          Plot.diagramVector(diagram2, device);
        else
          Plot.diagram(diagram2[1], device);
        end if;

    end bode;

    end Plot;

    encapsulated package Analysis

    end Analysis;

    end DiscreteTransferFunction;