Example 1

Description of the problem

Consider the synthetic continuous bilinear dynamical system given in its state-space form by

\[\begin{split}\dot{\boldsymbol{x}}(t) &= A_c\boldsymbol{x}(t) + \sum_{i=1}^2N_{c_i}\boldsymbol{x}(t)u_i(t) + B_c\boldsymbol{u}(t)\\ \boldsymbol{y}(t) &= C\boldsymbol{x}(t) + D\boldsymbol{u}(t)\end{split}\]

with

\[\begin{split}A_c = \begin{bmatrix} -1 & 0\\ 1 & -2 \end{bmatrix}, \quad N_{c_1} = \begin{bmatrix} 0 & 0\\ 1 & 1 \end{bmatrix}, \quad N_{c_2} = \begin{bmatrix} 1 & 1\\ 0 & 0 \end{bmatrix}, \quad B_{c} = \begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix}, \quad C = \begin{bmatrix} 0 & 1 \end{bmatrix}, \quad D = \begin{bmatrix} 0 & 0 \end{bmatrix}.\end{split}\]

The procedure in Python using the package systemID is highlighted below, with \(N_1 = N_2 = 10\), a frequency of acquisition \(f = 20\) Hz for a total time of \(5\) seconds. As is the case for linear systems, the realized system matrices are not unique, because the state space description is not unique. However, the input/output mapping should be unique and the linear part of the identified system matrix should have the same eigenvalues as the true system matrix. The errors in the system matrix eigenvalues (between true and identified) are

\[\begin{split}\left|\left| \lambda\left(A_c\right) - \lambda\left(\hat{A}_c\right)\right|\right| & \simeq 10^{-13}\\ \left|\left| \lambda\left(N_{c_1}\right) - \lambda\left(\hat{N}_{c_1}\right)\right|\right| & \simeq 10^{-12}\\ \left|\left| \lambda\left(N_{c_2}\right) - \lambda\left(\hat{N}_{c_2}\right)\right|\right| & \simeq 10^{-12}\end{split}\]

The identified system was subject to some test inputs and the response from the true system to the same test inputs was performed. The test inputs applied to the plants are

\[\begin{split}\boldsymbol{u}(t) = \begin{bmatrix} \sin(7t)\\ \cos(10t) \end{bmatrix}.\end{split}\]

Code using systemID

## Imports
import numpy as np
from scipy import interpolate
from systemID.ClassesGeneral.ClassSystem import ContinuousBilinearSystem
from systemID.ClassesGeneral.ClassSignal import ContinuousSignal, OutputSignal, DiscreteSignal, subtract2Signals
from systemID.ClassesGeneral.ClassExperiments import Experiments
from systemID.ClassesSystemID.ClassBilinear import BilinearSystemID


## Define continuous time system matrices and parameters
state_dimension = 2
input_dimension = 2
output_dimension = 1

def A(t):
    return np.array([[-1, 0], [1, -2]])
def N(t):
    return np.array([[0, 0, 1, 1], [1, 1, 0, 0]])
def B(t):
    return np.array([[1, 0], [0, 1]])
def C(t):
    return np.array([[0, 1]])
def D(t):
    return np.zeros([output_dimension, input_dimension])


## Identification parameters
p = 10


## Signal parameters
frequency = 20
dt = 1 /frequency
total_time_training = 6
total_time_testing = 5
number_steps_training = total_time_training * frequency + 1
number_steps_testing = total_time_testing * frequency + 1
tspan_training = np.linspace(0, total_time_training, number_steps_training)
tspan_testing = np.linspace(0, total_time_testing, number_steps_testing)


## Create system
x0 = np.zeros(state_dimension)
system = ContinuousBilinearSystem(state_dimension, input_dimension, output_dimension, [(x0, 0)], 'Nominal system', A, N, B, C, D)


## Test signal
def u(t):
    return np.array([np.sin(7*t), np.cos(10*t)])
test_signal = ContinuousSignal(input_dimension, signal_shape='External', u=u)
test_signal_d = DiscreteSignal(input_dimension, total_time_testing, frequency, signal_shape='External', data=u(tspan_testing))
true_output = OutputSignal(test_signal, system, tspan=tspan_testing)


## Experiments
N1 = 10
data_inputs_1 = np.zeros([input_dimension, number_steps_training, N1])
data_inputs_1[:, 0, :] = np.random.randn(input_dimension, N1)

N2 = 10
data_inputs_2 = np.zeros([input_dimension, number_steps_training, N2, N2])
for i in range(N2):
    random_input = data_inputs_1[:, 0, i]
    data_inputs_2[:, 0, :, i] = data_inputs_1[:, 0, 0:N2]
    data_inputs_2[:, 1, :, i] = np.outer(random_input, np.ones(N2))

inputs_1 = []
inputs_2 = []
systems = []

for i in range(N1):
    inputs_1.append(ContinuousSignal(input_dimension, signal_shape='External', u=interpolate.interp1d(tspan_training, data_inputs_1[:, :, i], kind='zero')))
    systems.append(system)

for i in range(N2):
    inputs_2.append([])
    for j in range(N2):
        inputs_2[-1].append(ContinuousSignal(input_dimension, signal_shape='External', u=interpolate.interp1d(tspan_training, data_inputs_2[:, :, j, i], kind='zero')))

experiments_1 = Experiments(systems, inputs_1, tspan=tspan_testing, total_time=total_time_testing, frequency=frequency)
experiments_2 = []
for i in range(N2):
    experiments_2.append(Experiments(systems[0:N2], inputs_2[i], tspan=tspan_testing, total_time=total_time_testing, frequency=frequency))


## Identification
bilinear = BilinearSystemID(experiments_1, experiments_2, state_dimension, dt, p=p)


## Identified system
x0_id = x0
identified_system = ContinuousBilinearSystem(state_dimension, input_dimension, output_dimension, [(x0_id, 0)], 'Identified system', bilinear.A, bilinear.N, bilinear.B, bilinear.C, bilinear.D)


## Test
identified_output = OutputSignal(test_signal, identified_system, tspan=tspan_testing)

Results

Output profiles obtained from the true and identified systems are compared below.

Alternative text

The singular value decomposition plot is displayed at different time instants.

Alternative text