Getting State Space Representation from system

[PauloMorgs]

Hello folks,

I would like to know how I would go on about getting the state space representation A,B,C and D of my model?

(Something like model.get_state_space_representation() from PyFMI)

Thanks!

[t-sommer]

This functionality does currently not exist in FMPy.

[t-sommer]

Can you give this code a try?

StateSpace.zip

model StateSpace

  parameter Real A[3,3] =
     [1.1, 2.1, 3.1;
      4.1, 5.1, 6.1;
      7.1, 8.1, 9.1];

  parameter Real B[3,3] =
     [1.2, 2.2, 3.2;
      4.2, 5.2, 6.2;
      7.2, 8.2, 9.2];

  parameter Real C[3,3] =
     [1.3, 2.3, 3.3;
      4.3, 5.3, 6.3;
      7.3, 8.3, 9.3];

  parameter Real D[3,3] =
     [1.4, 2.4, 3.4;
      4.4, 5.4, 6.4;
      7.4, 8.4, 9.4];

  input Real u[3];

  output Real y[3];

  Real x[3](start={0, 0, 0}, fixed=true);

equation 

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

end StateSpace;

from pathlib import Path

from tempfile import TemporaryDirectory

from _ctypes import byref
from numpy.typing import NDArray

from fmpy import read_model_description, instantiate_fmu, extract
from fmpy.fmi2 import FMU2Model, fmi2ValueReference, fmi2Real, fmi2EventInfo, _FMU2
import numpy as np

from fmpy.model_description import ModelDescription


def linearize(fmu: _FMU2, model_description: ModelDescription) -> tuple[NDArray, NDArray, NDArray, NDArray]:
    """
    Computes the linear state-space matrices (A, B, C, D) for an FMU using directional derivatives.

    Args:
        fmu: The FMU instance to linearize.
        model_description: The model description containing variables and structure.

    Returns:
        The state-space matrices:
            - A: State matrix (n x n)
            - B: Input matrix (n x m)
            - C: Output matrix (r x n)
            - D: Feedthrough matrix (r x m)
        where n is the number of states, m is the number of inputs, and r is the number of outputs.
    """

    vr_dx: list[int] = []
    vr_x: list[int] = []

    for unknown in model_description.derivatives:
        vr_dx.append(unknown.variable.valueReference)
        vr_x.append(unknown.variable.derivative.valueReference)

    vr_y: list[int] = []

    for unknown in model_description.outputs:
        vr_y.append(unknown.variable.valueReference)

    vr_u: list[int] = []

    for variable in model_description.modelVariables:
        if variable.causality == "input" and variable.variability == "continuous":
            vr_u.append(variable.valueReference)

    n = len(vr_x)  # number of continuous states
    m = len(vr_u)  # number of inputs
    r = len(vr_y)  # number of outputs

    A = np.zeros(shape=(n, n))
    B = np.zeros(shape=(n, m))
    C = np.zeros(shape=(r, n))
    D = np.zeros(shape=(r, m))

    vr_known = (fmi2ValueReference * 1)()
    seed = (fmi2Real * 1)(1)

    # A
    vr_unknown = (fmi2ValueReference * n)(*vr_dx)
    d_unknown = (fmi2Real * n)()

    for i, vr in enumerate(vr_x):
        vr_known[0] = vr
        fmu.fmi2GetDirectionalDerivative(
            fmu.component,
            vr_unknown,
            len(vr_unknown),
            vr_known,
            len(vr_known),
            seed,
            d_unknown,
        )
        A[:, i] = d_unknown

    # B
    vr_unknown = (fmi2ValueReference * n)(*vr_dx)

    for i, vr in enumerate(vr_u):
        vr_known[0] = vr
        fmu.fmi2GetDirectionalDerivative(
            fmu.component,
            vr_unknown,
            len(vr_unknown),
            vr_known,
            len(vr_known),
            seed,
            d_unknown,
        )
        B[:, i] = d_unknown

    # C
    vr_unknown = (fmi2ValueReference * r)(*vr_y)
    d_unknown = (fmi2Real * r)()

    for i, vr in enumerate(vr_x):
        vr_known[0] = vr
        fmu.fmi2GetDirectionalDerivative(
            fmu.component,
            vr_unknown,
            len(vr_unknown),
            vr_known,
            len(vr_known),
            seed,
            d_unknown,
        )
        C[:, i] = d_unknown

    # D
    vr_unknown = (fmi2ValueReference * r)(*vr_y)

    for i, vr in enumerate(vr_u):
        vr_known[0] = vr
        fmu.fmi2GetDirectionalDerivative(
            fmu.component,
            vr_unknown,
            len(vr_unknown),
            vr_known,
            len(vr_known),
            seed,
            d_unknown,
        )
        D[:, i] = d_unknown

    return A, B, C, D

tempdir =  TemporaryDirectory()

unzipdir = extract("StateSpace.fmu", tempdir.name)

model_description = read_model_description(unzipdir)

vr = dict((variable.name, variable.valueReference) for variable in model_description.modelVariables)

fmu: FMU2Model = instantiate_fmu(
    unzipdir=unzipdir,
    model_description=model_description,
    fmi_type="ModelExchange",
    fmi_call_logger=print,
)

fmu.fmi2SetupExperiment(fmu.component, False, 0.0, 0.0, True, 1.0)
fmu.fmi2EnterInitializationMode(fmu.component)
fmu.fmi2ExitInitializationMode(fmu.component)
eventInfo = fmi2EventInfo()
fmu.fmi2NewDiscreteStates(fmu.component, byref(eventInfo))
fmu.fmi2EnterContinuousTimeMode(fmu.component)

A, B, C, D = linearize(fmu, model_description)

fmu.freeInstance()

print(f"{A=}")
print(f"{B=}")
print(f"{C=}")
print(f"{D=}")

[t-sommer]

@casella, you can convert the ABCD matrices to sparse representations using scipy.sparse.csr_array and scipy.sparse.bsr_array.

[casella]

Thank @t-sommer, we'll try this out and report next week.