# Low Parametric Closed-Loop Sensitivity Realizations using Fixed-Point and Floating-Point Arithmetic

Thibault Hilaire, Philippe Chevrel, James F Whidborne

► To cite this version:

Thibault Hilaire, Philippe Chevrel, James F Whidborne. Low Parametric Closed-Loop Sensitivity Realizations using Fixed-Point and Floating-Point Arithmetic. European Control Conference, Jul 2007, Kos, Greece. pp.6, 10.23919/ECC.2007.7068839. [https://imt-atlantique.hal.science/hal-01317344](https://imt-atlantique.hal.science/hal-01317344)

HAL Id: hal-01317344

https://imt-atlantique.hal.science/hal-01317344

Submitted on 8 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Low Parametric Closed-Loop Sensitivity Realizations using Fixed-Point and Floating-Point Arithmetic

T. Hilaire*, P. Chevrel†, and J.F. Whidborne‡,

*IRISA UMR CNRS 6074, 6 rue Kérampont

22305 LANNION, FRANCE

†IRCCyN, UMR CNRS 6597, 1 rue de la Noë

44321 NANTES, FRANCE

‡Department of Aerospace Science, Cranfield University

Bedfordshire MK43 0AL, U.K.

**Abstract**— The Specialized Implicit Form provides a general framework for the analysis and design of digital controller implementations with minimal finite wordlength effects. This paper proposes a measure of the closed-loop transfer function sensitivity to finite wordlength effects that is generalized for both fixed-point and floating-point arithmetic and can be used with the Specialized Implicit Form to analyze the effects of quantization and rounding on the parameters of a digital controller implementation. The measure is computationally tractable and hence amenable to solving the problem of minimizing the parametric sensitivity FWL effect. Furthermore, the sensitivity to the rounding of each individual parameter can be easily obtained. The use of the measure is illustrated with examples.

**Index Terms**— Finite-wordlength effects, Implementation, Fixed-point arithmetic, Floating-point arithmetic, Controller realization

## I. INTRODUCTION

When control systems are implemented in digital hardware, rounding and quantization occurs on the variables and constants in the controller resulting from the finite-precision nature of the number storage in the computing device. This arises because a number must be represented with a finite wordlength (FWL). There are two main effects of this finite-precision (often known as the *FWL effects*). The first is the addition of noise into the system resulting from the rounding of variables before and after each arithmetic operation - the “round-off noise”. The second is the degradation in the performance and/or the stability resulting from quantization of the controller coefficients/parameters, known as “parametric sensitivity”, about which this paper is concerned. For most low-order controllers, the FWL effects are insignificant, but for higher-order controllers, particularly when fast sampling is used, the FWL effects can become significant. However, it is well-known that the FWL effects are dependent upon the controller realization, hence there has been a great deal of work in determining realizations that minimize the FWL effects in some sense, e.g. [1], [2], [3]. It is also well-known that the FWL effects are dependent on the choice

of dynamic operator. The $\delta$-operator, for example, generally has much better numerical properties than the usual delay operator, $q$, for control systems with fast sampling, e.g. [4]. Unsurprisingly, the FWL effects are also dependent upon the choice of arithmetic, either fixed-point or floating-point.

The problem of addressing the optimal realization for minimal FWL effects is usually addressed in the state space, e.g. [5], [1], [3]. Briefly, if the controller is $K(s) = C(\sigma I - A)^{-1}B + D$ where $\sigma$ is usually the transform of the $\delta$ or $q$-operator, the problem is to search over the set

$$ \{ CT(\sigma I - T^{-1}AT)^{-1}TB + D : T \text{ a non-singular matrix} \} $$

to find a $T$ and hence controller realization that is insensitive to FWL effects. The limitation of this approach is that (i) there are many realizations that cannot be expressed in standard state space form and (ii) the search is confined to a single operator. The $\delta$-operator is more complex to implement than the $q$-operator, so in some circumstances, it may be best to have a mix of operators. These limitations may be overcome by use of the *Specialized Implicit Form* for the controller [6], [7], [8].

This paper proposes a measure of the sensitivity of the closed-loop transfer function to the controller parameter rounding that is based on a measure proposed by [1] and is suitable for use with the specialized implicit form. The measure is generalized for either fixed or floating-point arithmetic and can be easily computed. In the next section, the specialized implicit form is introduced. In section III, the quantization of the controller coefficients is discussed, and in the following section, the sensitivity measure is proposed and the means of its evaluation provided. In section V, the problem of determining a realization with minimal parametric sensitivity is posed and some examples provided.

## II. A UNIFYING FRAMEWORK

### A. The *Specialized Implicit Form*

Many controller forms, such as lattice filters and $\delta$-operator controllers, use intermediate variables, and hence cannot be expressed in the traditional state-space form. The Specialized Implicit Form allows a description, in a

Email: thibault.hilaire@irisa.fr, philippe.chevrel@emn.fr and j.f.whidborne@cranfield.ac.uk
single equation, of almost any implementation. Furthermore, it provides an explicit description of the parametrization, and allows the analysis of the FWL effects, but is still a macroscopic description. The description takes the form of an implicit state-space system [9] and is given by

$$ \begin{pmatrix} J & 0 & 0 \\ -K & I_n & 0 \\ -L & 0 & I_p \end{pmatrix} \begin{pmatrix} T(k+1) \\ X(k+1) \\ Y(k) \end{pmatrix} = \begin{pmatrix} 0 & M & N \\ 0 & P & Q \\ 0 & R & S \end{pmatrix} \begin{pmatrix} T(k) \\ X(k) \\ U(k) \end{pmatrix} \quad (1) $$

where

* $J \in \mathbb{R}^{l \times l}$, $K \in \mathbb{R}^{n \times l}$, $L \in \mathbb{R}^{p \times l}$, $M \in \mathbb{R}^{l \times n}$, $N \in \mathbb{R}^{l \times m}$, $P \in \mathbb{R}^{n \times n}$, $Q \in \mathbb{R}^{n \times m}$, $R \in \mathbb{R}^{p \times n}$, $S \in \mathbb{R}^{p \times m}$, $T(k) \in \mathbb{R}^l$, $X(k) \in \mathbb{R}^n$, $U(k) \in \mathbb{R}^m$ and $Y(k) \in \mathbb{R}^p$,

* matrix $J$ is lower triangular with 1's on the main diagonal.

* $T(k+1)$ is the intermediate variable in the calculations of step $k$ (the column of 0's in the second matrix shows that $T(k)$ is not used for the calculation at step $k$ - this characterizes the concept of an intermediate variable),

* $X(k+1)$ is the stored state-vector ($X(k)$ is effectively stored from one step to the next, in order to compute $X(k+1)$ at step $k$).

$T(k+1)$ and $X(k+1)$ form the descriptor-vector: $X(k+1)$ is stored from one step to the next, while $T(k+1)$ is computed and used within one time step.

It is implicitly assumed throughout the paper that the computations associated with the realization (1) are executed in row order giving the following algorithm:

$$ [i] \quad JT(k+1) \leftarrow MX(k) + NU(k) $$

$$ [ii] \quad X(k+1) \leftarrow KT(k+1) + PX(k) + QU(k) $$

$$ [iii] \quad Y(k) \leftarrow LT(k+1) + RX(k) + SU(k) $$

Note that in practice, steps [ii] and [iii] could be exchanged to reduce the computational delay. Also note that because the computations are executed in row order and $J$ is lower triangular with 1's on the main diagonal, there is no need to compute $J^{-1}$.

Equation (1) is equivalent in infinite precision to the classical state-space form

$$ \begin{pmatrix} T(k+1) \\ X(k+1) \\ Y(k) \end{pmatrix} = \begin{pmatrix} 0 & J^{-1}M & J^{-1}N \\ 0 & A_Z & B_Z \\ 0 & C_Z & D_Z \end{pmatrix} \begin{pmatrix} T(k) \\ X(k) \\ U(k) \end{pmatrix} \quad (2) $$

with $A_Z \in \mathbb{R}^{n \times n}$, $B_Z \in \mathbb{R}^{n \times m}$, $C_Z \in \mathbb{R}^{p \times n}$ and $D_Z \in \mathbb{R}^{p \times m}$ where

$$ A_Z = KJ^{-1}M + P, \quad B_Z = KJ^{-1}N + Q, \quad (3) $$

$$ C_Z = LJ^{-1}M + R, \quad D_Z = LJ^{-1}N + S. \quad (4) $$

Note that (2) corresponds to a different parametrization than (1). The system transfer function is given by

$$ H(z) = C_Z(zI_n - A_Z)^{-1}B_Z + D_Z. \quad (5) $$

## B. Definitions

To complete the framework, the following definitions are required. For further details, see [8], [10].

**Definition 1** A realization, $\mathcal{R}$, is defined by the specific set of matrices $J, K, L, M, N, P, Q, R$ and $S$ used to describe a realization with the implicit form of (1):

$$ \mathcal{R} \triangleq (J, K, L, M, N, P, Q, R, S). \quad (6) $$

**Remark 1** $\mathcal{R}$ can also be defined by the matrix $Z \in \mathbb{R}^{(l+n+p) \times (l+n+m)}$

$$ Z \triangleq \begin{pmatrix} -J & M & N \\ K & P & Q \\ L & R & S \end{pmatrix} \quad (7) $$

and the dimensions $l, m, n$ and $p$, so $\mathcal{R}$ could be defined by $\mathcal{R} := (Z, l, m, n, p)$.

**Definition 2** $\mathcal{R}_H$ denotes the set of realizations with transfer function $H$. These realizations are said to be equivalent.

In order to encompass realizations with some special structure (q-operator state-space, $\delta$-operator state-space, direct form, cascade, lattice filters, etc.), we define a set of realizations that possess a particular structure.

**Definition 3** A structuration $\mathcal{S}$ is a set of realizations having a common structure: some coefficients or some dimensions are fixed a priori.

Some examples of common structurations are given in the next section.

**Definition 4** $\mathcal{R}_H^\mathcal{S}$ is the set of equivalent structured realizations. Realizations from $\mathcal{R}_H^\mathcal{S}$ are structured according to $\mathcal{S}$ and have a transfer function $H$. Hence

$$ \mathcal{R}_H^\mathcal{S} \triangleq \mathcal{R}_H \cap \mathcal{S}. \quad (8) $$

Some coefficients of $Z$ are, for example, always set to unity or zero, and are hence not subject to rounding. These are not 'significant coefficients' and hence are not included in the parametrization.

**Definition 5** A parametrization of a realization $\mathcal{R}$ is the set of coefficients of $Z$ that are significant for the realization.

### C. Some examples

The following specialized implicit form describes the realization of a state-space controller using the $\delta$-operator

$$ \begin{pmatrix} I_n & 0 & 0 \\ -\Delta I_n & I_n & 0 \\ 0 & 0 & I_p \end{pmatrix} \begin{pmatrix} T(k+1) \\ X(k+1) \\ Y(k) \end{pmatrix} = \begin{pmatrix} 0 & A_\delta & B_\delta \\ 0 & I_n & 0 \\ 0 & C_\delta & D_\delta \end{pmatrix} \begin{pmatrix} T(k) \\ X(k) \\ U(k) \end{pmatrix} $$

So, the $\delta$-structuration $\mathcal{S}_\delta$ is formally defined by

$$ \mathcal{S}_\delta = \left\{ \begin{array}{l} \mathcal{R} := (I_n, \Delta I_n, 0, A_\delta, B_\delta, I_n, 0, C_\delta, D_\delta) \\ \quad \forall m \in \mathbb{N}, n \in \mathbb{N}, p \in \mathbb{N} \\ \quad \forall \Delta \in \mathbb{R}^+, A_\delta \in \mathbb{R}^{n \times n}, B_\delta \in \mathbb{R}^{n \times m} \\ \quad \forall C_\delta \in \mathbb{R}^{p \times n}, D_\delta \in \mathbb{R}^{p \times m} \end{array} \right\} \quad (9) $$
The cascade form is a common realization for filter implementation. It generally has good FWL properties compared to the direct forms. For cascade form, the filter is decomposed into a number of lower order (usually first and second-order) transfer function blocks connected in series. For the next example, we consider two standard q-operator filter blocks connected in series as shown in fig. 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Cascade form</figcaption></figure>

If the two state-space realizations $\mathcal{R}_1$ and $\mathcal{R}_2$ are defined by $(A_1, B_1, C_1, D_1)$ and $(A_2, B_2, C_2, D_2)$, then cascading $\mathcal{R}_1$ with $\mathcal{R}_2$ leads to the following realization

$$Z = \begin{pmatrix} -I & C_1 & 0 & D_1 \\ 0 & A_1 & 0 & B_1 \\ B_2 & 0 & A_2 & 0 \\ D_2 & 0 & C_2 & 0 \end{pmatrix} \quad (10)$$

from which definition of the structuration $\mathcal{S}$ immediately follows. The output of $\mathcal{R}_1$ is computed in the intermediate variable, and used as the input of $\mathcal{R}_2$.

The main point is that if we consider the equivalent state-space realization, with parameters

$$A = \begin{pmatrix} A_1 & 0 \\ B_2 C_1 & A_2 \end{pmatrix}, \quad B = \begin{pmatrix} B_1 \\ B_2 D_1 \end{pmatrix}, \quad (11) \\ C = \begin{pmatrix} D_2 C_1 & C_2 \end{pmatrix}, \quad D = D_2 D_1,$$

the parametrization is not the one used in the computations.

For a given form, it is generally straightforward to define the structuration. A number of other examples are given in [8] and [6].

### III. COEFFICIENT QUANTIZATION

A coefficient's quantization depends both on their value and their representation.

Firstly if the value of a coefficient is such that it will be quantized without error, then that parameter makes no contribution to the overall coefficient sensitivity. Hence we introduce weighting matrices $W_J$ to $W_S$ (and also $W_Z$) respectively associated with matrices $J$ to $S$ of a realization, such that

$$(W_X)_{i,j} \triangleq \begin{cases} 0 & \text{if } X_{i,j} \text{ is exactly implemented,} \\ 1 & \text{otherwise.} \end{cases} \quad (12)$$

Secondly, different representation schemes may be considered. Here we consider both fixed-point and floating-point representations of coefficients expressed using $\beta$ bits.

A fixed-point coefficient $x$ is represented by $(-1)^s \cdot N \cdot 2^{-\beta_f}$, where $s \in \{0, 1\}$, $N$ is an integer coded with $\beta_g$ bits and $\beta_f$ an integer (not stored in the representation) such that $\beta_g + \beta_f + 1 = \beta$. The quantized $x^\dagger$ of $x$ is such that

$$|x - x^\dagger| < 2^{-(\beta_f + 1)}. \quad (13)$$

A floating-point coefficient is represented by $(-1)^s \cdot N \cdot 2^e$ where $s \in \{0, 1\}$ (or $w \in [0.5, 1]$ for a normalized floating-point representation) and $e$ is an integer coded with $\beta_e$ bits<fnref n="1" /> ($\beta_e + \beta_w + 1 = \beta$). The quantized $x^\dagger$ of $x$ is, in this case, such that

$$|x - x^\dagger| < x \cdot 2^{-(\beta_w + 1)}. \quad (14)$$

The choice of $\beta_f$ and $e$ can be unique for each coefficient ($e = \lceil \log_2 |x| \rceil$ and $\beta_f = \beta - 1 - \lceil \log_2 |x| \rceil$, where $\lceil \cdot \rceil$ is the ceiling operator). Alternatively, $\beta_f$ and $e$ are defined for a group of coefficients (in order to reduce the required bit-shifts and the subsequent computational cost). This defines the block-fixed-point and block-floating-point schemes. Following [11], we introduce the generalized dynamic range bit $\beta_r$ ($\beta_r = \beta_g$ or $\beta_e$) and the precision bit length $\beta_p$ ($\beta_p = \beta_f$ or $\beta_w$).

Usually, the blocks used in block-representation correspond to the matrices $J$ to $S$, but there is no necessity for this, and blocks can be chosen at will. To define the blocks of a realization $\mathcal{R} := (Z, l, m, n, p)$, we introduce the matrix $\eta_Z$ such that

$$(\eta_Z)_{i,j} \triangleq \begin{cases} \text{the largest absolute value of} \\ \text{the block in which } Z_{i,j} \text{ resides.} \end{cases} \quad (15)$$

This allows a completely general definition of the blocks. Thus there could be just a single unique block, or every block could consist of only one coefficient. For example, denoting $\mathcal{E}_{a,b} \in \mathbb{R}^{a \times b}$ as a matrix of 1s and

$$\|X\|_{\max} \triangleq \max_{i,j} |X_{i,j}|, \quad (16)$$

then using a block-representation corresponding to the matrices $J$ to $S$ gives

$$\eta_Z = \begin{pmatrix} \|J\|_{\max} & \mathcal{E}_{l,l} & \|M\|_{\max} & \mathcal{E}_{l,n} & \|N\|_{\max} & \mathcal{E}_{l,m} \\ \|K\|_{\max} & \mathcal{E}_{n,l} & \|P\|_{\max} & \mathcal{E}_{n,n} & \|Q\|_{\max} & \mathcal{E}_{n,m} \\ \|L\|_{\max} & \mathcal{E}_{p,l} & \|R\|_{\max} & \mathcal{E}_{p,n} & \|S\|_{\max} & \mathcal{E}_{p,m} \end{pmatrix}.$$

With a single unique block for $Z$ we get

$$\eta_Z = \|Z\|_{\max} \mathcal{E}_{l+n+p,l+n+m}, \quad (17)$$

and for one block per coefficient we get

$$(\eta_Z)_{i,j} = |Z_{i,j}|. \quad (18)$$

**Proposition 1** During the quantization process, $Z$ is perturbed to $Z + r_Z \times \Delta$ where

$$r_Z \triangleq \begin{cases} W_Z & \text{for fixed-point representation,} \\ 2\eta_Z \times W_Z & \text{for floating-point representation,} \end{cases} \quad (19)$$

$\Delta$ is a matrix dependant on the $\beta_p$ precision bit length, and $\times$ denotes the Schur product. If $\beta_{p,i,j}$ is the precision bit-length of $Z_{i,j}$, then

$$|\Delta_{i,j}| < 2^{-(\beta_{p,i,j} + 1)}. \quad (20)$$

<footnote n="1">¹The difference with fixed-point is that $e$ is coded with $\beta_e$ bits and can changed. With fixed-point, $\beta_f$ is fixed and implicit.</footnote>
**Proof:** The proof comes from the quantized error expressed in (13) and (14). ■

**Remark 2** With this formalism for the different representation schemes, note that the choice of the scale parameter ($e$ or $\beta_f$) is defined for each coefficient ($e_{i,j} = \lceil \log_2 |\eta_{Z_{i,j}}| \rceil$ and $\beta_{f_{i,j}} = \beta - 1 - \lceil \log_2 |\eta_{Z_{i,j}}| \rceil$) and that it is also possible to define the minimum bit length $\beta$ to code each coefficient without overflow or underflow [8].

#### IV. CLOSED-LOOP SENSITIVITY MEASURES

In order to determine the optimal realization, some measure of the closed-loop sensitivity to parametric perturbations is required. A fair number of these have been proposed over the years. Ideally, the chosen measure should be computationally tractable but reasonably representative of the actual perturbations that occur in implementation. Probabilistic measures have been proposed, e.g. [12], as well as small-gain measures [13]. A number of measures of the closed-loop pole sensitivity have been proposed [14], [3], [15], [16]. Here though, we use a closed-loop transfer function sensitivity measure which extends that originally proposed by [1].

##### A. Problem statement

<figure><img src="image_3.png" /><figcaption>Fig. 2. Closed-loop system considered</figcaption></figure>

Consider a plant, $\mathcal{P}$, with controller, $\mathcal{C}$, in the standard form shown in fig. 2, where $W(k) \in \mathbb{R}^{p_1}$ and $Z(k) \in \mathbb{R}^{m_1}$ are the exogenous input and output respectively and $U(k) \in \mathbb{R}^{p_2}$ and $Y(k) \in \mathbb{R}^{m_2}$ are the plant control and measured signals respectively.

The plant $\mathcal{P}$ is defined by

$$ \mathcal{P} := \left( \begin{array}{c|cc} A & B_1 & B_2 \\ \hline C_1 & D_{11} & D_{12} \\ C_2 & D_{21} & 0 \end{array} \right) \qquad (21) $$

where $A \in \mathbb{R}^{n_p \times n_p}$, $B_1 \in \mathbb{R}^{n_p \times p_1}$, $B_2 \in \mathbb{R}^{n_p \times p_2}$, $C_1 \in \mathbb{R}^{m_1 \times n_p}$, $C_2 \in \mathbb{R}^{m_2 \times n_p}$, $D_{11} \in \mathbb{R}^{m_1 \times p_1}$, $D_{12} \in \mathbb{R}^{m_1 \times p_2}$, $D_{21} \in \mathbb{R}^{m_2 \times p_1}$ and $D_{22} \in \mathbb{R}^{m_2 \times p_2}$ is assumed to be zero only in order to simplify the expressions. The controller is defined by a realization $\mathcal{C} := (Z, l, m, n, p_2)$ with transfer function $C_Z(zI - A_Z)^{-1}B_Z + D_Z$.

The closed-loop system $\mathcal{S}$ is then given by

$$ \mathcal{S} = Fl(\mathcal{P}, \mathcal{C}) := \left( \begin{array}{c|c} \bar{A} & \bar{B} \\ \hline \bar{C} & \bar{D} \end{array} \right) \qquad (22) $$

where $Fl(\cdot, \cdot)$ is the well-known lower fractional transform operation [17] and where $\bar{A} \in \mathbb{R}^{n_p+n \times n_p+n}$, $\bar{B} \in \mathbb{R}^{n_p+n \times p_1}$, $\bar{C} \in \mathbb{R}^{m_1 \times n_p+n}$ and $\bar{D} \in \mathbb{R}^{m_1 \times p_1}$ such that

$$ \bar{A} = \begin{pmatrix} A + B_2 D_Z C_2 & B_2 C_Z \\ B_Z C_2 & A_Z \end{pmatrix} \qquad (23) $$

$$ \bar{B} = \begin{pmatrix} B_1 + B_2 D_Z D_{21} \\ B_Z D_{21} \end{pmatrix} \qquad (24) $$

$$ \bar{C} = (C_1 + D_{12} D_Z C_2 \quad D_{12} C_Z) \qquad (25) $$

$$ \bar{D} = D_{11} + D_{12} D_Z D_{21} \qquad (26) $$

The closed-loop transfer function is

$$ \bar{H} : z \to \bar{C} (zI - \bar{A})^{-1} \bar{B} + \bar{D} \qquad (27) $$

##### B. Transfer function sensitivity

In order to evaluate how much the digital approximation of the controller's coefficients (due to FWL implementation) affects the closed-loop transfer function, the transfer function sensitivity $\frac{\partial \bar{H}}{\partial Z}$ can be used.

Let $\bar{H}^{\dagger} \triangleq \bar{H}|_{Z+r_Z \times \Delta}$ denote the closed-loop transfer function $\bar{H}$ perturbed by the quantization process ($Z$ is perturbed in $Z+r_Z \times \Delta$ according to proposition 1). Then, for the Single Input Single Output (SISO) case, $\forall z \in \mathbb{C}$

$$ \bar{H}^{\dagger}(z) - \bar{H}(z) = \sum_{i,j} \Delta_{i,j} \left. \frac{\partial \bar{H}^{\dagger}(z)}{\partial \Delta} \right|_{\Delta=0} + o\left(\|\Delta\|_{\max}\right) \qquad (28) $$

and

$$ \| \bar{H}^{\dagger} - \bar{H} \|_{2} \lesssim \| \Delta \|_{\max} \left\| \left. \frac{\partial \bar{H}^{\dagger}}{\partial \Delta} \right|_{\Delta=0} \right\|_{2} \qquad (29) $$

where $\|\cdot\|_2$ is the $L_2$-norm. It is easy to see that

$$ \left. \frac{\partial \bar{H}^{\dagger}}{\partial \Delta} \right|_{\Delta=0} = \frac{\partial \bar{H}}{\partial Z} \times r_Z \qquad (30) $$

and so (29) leads to the following transfer function sensitivity measure:

**Definition 6** Let consider a realization $\mathcal{R} := (Z, l, m, n, p)$ with representation matrix $r_Z$. The closed-loop transfer function sensitivity, with respect to all the non-trivial coefficients of $\mathcal{R}$ is defined in the SISO case by

$$ \bar{M}_{L_2}^W \triangleq \left\| \frac{\partial \bar{H}}{\partial Z} \times r_Z \right\|_2^2 \qquad (31) $$

**Remark 3** This measure must be linked to the open-loop transfer function sensitivity $\left\| \frac{\partial \bar{H}}{\partial Z} \times r_Z \right\|_2^2$ previously defined in [7], [18], [10] that derives from Gevers' definition [1].

**Remark 4** From (31) and (29), it is possible to ensure that the closed-loop transfer function perturbation is smaller than a certain constant in an $L_2$-norm sense. It is also possible to include a frequency weighting to emphasize certain frequency ranges. Furthermore, it is also possible to use an $H_{\infty}$-norm to ensure that the closed-loop degradation is constrained over a given frequency range.
This measure can be extended to the Multiple Input Multiple Output (MIMO) case. However, it is also useful to be able to consider the contribution of each coefficient to the overall sensitivity. The closed-loop transfer function sensitivity matrix, denoted by $\frac{\delta \bar{H}}{\delta \bar{Z}}$, is the matrix of the $L_2$-norm of the sensitivity of the transfer function $\bar{H}$ with respect to each coefficient $Z_{i,j}$. It is defined by

$$ \left( \frac{\delta \bar{H}}{\delta \bar{Z}} \right)_{i,j} \triangleq \left\| \frac{\partial \bar{H}}{\partial Z_{i,j}} \right\|_2 \qquad (32) $$

and allows the evaluation of the overall impact of each coefficient. It can be used to evaluate the overall sensitivity. From the properties of $L_2$-norms, we have

$$ \left\| \frac{\delta \bar{H}}{\delta \bar{Z}} \right\|_F = \left\| \frac{\partial \bar{H}}{\partial \bar{Z}} \right\|_2 \qquad (33) $$

where $\|\cdot\|_F$ is the Frobenius norm.

Definition 6 can now be extended to the general case (MIMO and SISO) :

**Definition 7** The closed-loop transfer function sensitivity measure is defined by

$$ \bar{M}_{L_2}^W \triangleq \left\| \frac{\delta \bar{H}}{\delta \bar{Z}} \times r_Z \right\|_F^2 \qquad (34) $$

The transfer function sensitivity $\frac{\partial \bar{H}}{\partial \bar{Z}}$ is given by the following proposition

**Proposition 2**

$$ \frac{\partial \bar{H}}{\partial \bar{Z}} = \left[ \begin{array}{c} \bar{C} (zI - \bar{A})^{-1} \bar{M}_1 + \bar{M}_2 \\ \bar{N}_1 (zI - \bar{A})^{-1} \bar{B} + \bar{N}_2 \end{array} \right] \circled{\text{}} \left[ \begin{array}{c} \bar{C} (zI - \bar{A})^{-1} \bar{M}_1 + \bar{M}_2 \\ \bar{N}_1 (zI - \bar{A})^{-1} \bar{B} + \bar{N}_2 \end{array} \right] \qquad (35) $$

where $\circled{\text{}}$ is an operator defined by

$$ A \circled{\text{}} B \triangleq \mathrm{Vec}(A) \cdot \left[ \mathrm{Vec}(B^\top) \right]^\top, \qquad (36) $$

$\mathrm{Vec}(\cdot)$ is the classical operator that vectorizes a matrix, and

$$ \bar{M}_1 = \begin{pmatrix} B_2 L J^{-1} & 0 & B_2 \\ K J^{-1} & I_n & 0 \end{pmatrix} \qquad (37) $$

$$ \bar{N}_1 = \begin{pmatrix} J^{-1} N C_2 & J^{-1} M \\ 0 & I_n \\ C_2 & 0 \end{pmatrix} \qquad (38) $$

$$ \bar{M}_2 = \begin{pmatrix} D_{12} L J^{-1} & 0 & D_{12} \end{pmatrix} \qquad (39) $$

$$ \bar{N}_2 = \begin{pmatrix} J^{-1} N D_{21} \\ 0 \\ D_{21} \end{pmatrix} \qquad (40) $$

*Proof:* The proof is based on the following lemma :

**Lemma 1** Let $X$ be a matrix in $\mathbb{R}^{p \times l}$ and $G, H$ be two transfer function with values respectively in $\mathbb{C}^{m \times p}$ and $\mathbb{C}^{l \times n}$. $G$ and $H$ are supposed to be independent w.r.t. $X$. Then

$$ \frac{\partial(GXH)}{\partial X} = G \circled{\text{}} H \qquad (41) $$

$$ \frac{\partial(GX^{-1}H)}{\partial X} = (GX^{-1}) \circled{\text{}} (X^{-1}H) \qquad (42) $$

The proofs can be found in [8]. ■

## V. OPTIMAL DESIGN

Since the closed-loop sensitivity measure depends on the realization chosen to numerically realize the controller, it is of interest to find, among the equivalent realizations, those with good closed-loop FWL properties.

### A. Equivalent realizations

In order to exploit the potential offered by the specialized implicit form in improving implementations, it is necessary to describe sets of equivalent system realizations. In [10], the Inclusion Principle introduced by Šiljak and Ikeda [19], [20] in the context of decentralized control, is extended to the Specialized Implicit Form in order to characterize equivalent classes of realizations. Although this extension gives the formal description of equivalent classes, it is of practical interest to consider only realizations with the same dimensions, where transformation from one realization to another is only a similarity transformation.

**Proposition 3** Consider a realization $\mathcal{R}_0 := (Z_0, l, m, n, p)$. All realizations $\mathcal{R} := (Z, l, m, n, p)$ such that

$$ Z = \begin{pmatrix} \mathcal{Y} & \\ & U^{-1} \end{pmatrix}_{I_p} Z_0 \begin{pmatrix} \mathcal{W} & \\ & U \end{pmatrix}_{I_m} \qquad (43) $$

with $U \in \mathbb{R}^{n \times n}$, $\mathcal{Y} \in \mathbb{R}^{l \times l}$ and $\mathcal{W} \in \mathbb{R}^{l \times l}$ non-singular, are equivalent.

**Remark 5** Given a realization $Z_0$ in the cascade form of (10), it is possible to characterize a subset of similarity transformation that preserves the cascade structure. The equivalent realizations with this particular structure are given by the particular similarity transform (specialization of eq. (43))

$$ Z = \begin{pmatrix} \mathcal{Y}^{-1} & & & \\ & U^{-1} & & \\ & & W^{-1} & \\ & & & I_p \end{pmatrix} Z_0 \begin{pmatrix} \mathcal{Y} & & & \\ & U & & \\ & & W & \\ & & & I_m \end{pmatrix} \qquad (44) $$

In the SISO case, the transformation matrix $\mathcal{Y}$ is only a scale factor between the first and the second stages of the realization, that can be considered to bound coefficients into acceptable limits. Note that in the MIMO case, there is an extra degree freedom.

### B. Optimal realizations problem

The problem of determining the realization that is best in some sense can be posed as follows:

**Problem 1 (Optimal realization problem)** Consider a transfer function $H$ and a sensitivity measure $\mathcal{J}$. The optimal design problem is to find the best realization $\mathcal{R}_{opt}$ with transfer function $H$ according to the criteria $\mathcal{J}$, that is

$$ \mathcal{R}_{opt} = \arg \min_{\mathcal{R} \in \mathcal{R}_H} \mathcal{J}(\mathcal{R}). \qquad (45) $$
Due to the size of $\mathcal{R}_H$, this problem generally cannot be solved practically. Hence the following problem is introduced to restrict the search to some particular structuration, each one is searched over using proposition 3.

### Problem 2 (Optimal structured realization problem)

Consider some structurations $(\mathcal{S}_i)_{1\le i\le N}$. The problem to find the optimal structured realization $\mathcal{R}_{opt}^\mathcal{S}$ is

$$ \mathcal{R}_{opt}^{\mathcal{S}} = \min_{1 \le i \le N} \left( \arg \min_{\mathcal{R} \in \mathcal{R}_H^{\mathcal{S}_i}} \mathcal{J}(\mathcal{R}) \right). \quad (46) $$

Since the measure $\mathcal{J}$ could be non-smooth and/or non-convex, the Adaptive Simulated Annealing (ASA) [21], [22] method has been chosen to solve Problem 2. This method has worked well for other optimal realization problems [15].

### C. Example

The first numerical example, used here to evaluate the closed-loop sensitivity under various parametrizations, is a SISO fluid power control system studied in [23], [13]. The discrete-time (sampled at 2 kHz) plant $\mathcal{P}$ is given by $(A_p, B_p, C_p)$ in (48), and transformed to the standard form.

The initial realization $\mathcal{R}_0 := (Z_0, 0, 4, 1, 1)$ of the controller $C$ is given in the controllable canonical form in equation (49). It is important to note that the coefficients are given with only 4 digits, but, due to the sensitivity of this example, this could be not sufficient to define correctly the system. Parameters that may be approximated by the quantization process required for implementation are shown in bold font, and the weighting matrix is built accordingly.

The following realizations are considered :

$Z_0$: Direct Form II : it corresponds to the canonical form

$Z_1$: optimal (according to $\bar{M}_{L_2}^W$) classical state-space realization

$Z_2$: Direct Form II with $\delta$-operator ($\Delta = 2^{-5}$)

$Z_3$: optimal $\delta$-realization

Each realization, and its sensitivity matrix $\frac{\delta\tilde{H}}{\delta Z}$, is presented in equations (49) to (52).

The values of the sensitivity measures for the various realizations are given in Table I. Only the fixed-point format, with block-representation corresponding to the matrices $J$ to $S$, is considered here.

The results obtained are consistent with existing results on open-loop sensitivity – it is interesting to note that realization $Z_1$ has a greater sensitivity than $Z_2$, which shows the extremely good numerical properties of the $\delta$-operator.

The second numerical example is an active control of longitudinal vehicle oscillations studied in [24]. One significant aspect of vehicle driveability is the attenuation of the first torsional mode (resonance in the elastic parts) which produces unpleasant (0 to 10 Hz) longitudinal oscillations of the vehicle, known as shuffle. They can be reduced by means of a controller acting on the engine torque.

The powertrain was modeled in continuous-time form, and a continuous-time $H_\infty$ optimal controller was designed [24].

<caption>TABLE I</caption>

<table><caption>EXAMPLE 1: CLOSED-LOOP SENSITIVITIES AND COMPUTATIONAL COST FOR DIFFERENT REALIZATIONS</caption><thead><tr><th>realization</th><th>M&#773;<sup>W</sup><sub>L<sub>2</sub></sub></th><th>Nb. operations</th></tr></thead><tbody><tr><td>Z<sub>0</sub></td><td>1.102e+17</td><td>8 + 9×</td></tr><tr><td>Z<sub>1</sub></td><td>1.334e+6</td><td>20 + 25×</td></tr><tr><td>Z<sub>2</sub></td><td>2.507e+4</td><td>12 + 13×</td></tr><tr><td>Z<sub>3</sub></td><td>33.765</td><td>24 + 29×</td></tr></tbody></table>

<caption>TABLE II</caption>

<table><caption>EXAMPLE 2: CLOSED-LOOP SENSITIVITIES AND COMPUTATIONAL COST</caption><thead><tr><th>realization</th><th>M&#773;<sup>W</sup><sub>L<sub>2</sub></sub></th><th>Nb. operations</th></tr></thead><tbody><tr><td>direct form II</td><td>6.802e+21</td><td>20 + 21×</td></tr><tr><td>M&#773;<sup>W</sup><sub>L<sub>2</sub></sub>-optimal state-space</td><td>6.328e+4</td><td>110 + 121×</td></tr><tr><td>optimal Observer-State-Feedback</td><td>8.441e+4</td><td>121 + 131×</td></tr></tbody></table>

The discretized model $P(z)$ is given by equations (53) and (54), and a discrete-time realization of the controller is given by (54) and (55) – it corresponds to an internally balanced realization.

The different forms studied here are :

* direct form II

* $\bar{M}_{L_2}^W$-optimal state-space

* Observer-State-Feedback

The Observer-State-Feedback form allows an enrichment of the observer model with a physical meaning but also because these states estimate the states of the physical system. Thus it improves the readability of the signals, and the initialization of the controller states is based on the physical states of the system. The Observer-State-Feedback form is given by

$$ \begin{cases} \hat{X}_{k+1} = A_p \hat{X}_k + B_p U_k + K_f (Y_k - C_p \hat{X}_k) \\ U_k = -K_c \hat{X}_k + Q (Y_k - C \hat{X}_k) \end{cases} \quad (47) $$

and can be written in the Specialized Implicit Form according to equation (56).

The transformation from the state-space form to the Observer-State-Feedback form requires the solution of a generalized Riccati equation [25]. The controller poles must be classified in three categories, which are the observation gain, the filter gain and the Youla parameter (static here). This repartition (according to some rules [18]) determines the parameters $K_f, K_c$ and $Q$. It is possible to numerically implement equation (47) in various ways, depending on the choice of the partition of the poles. Here, with 20 poles, there exist 184756 partitions, but only 120 are actually possible. The optimal realization is chosen from these realizations.

Table II shows the different sensitivity values. Note that the Observer-State-Feedback form is not significantly more sensitive than the state-space form, and with only a small increase in the number of arithmetic operations.

From Tables I and II it is clear that the optimal realizations, which are fully parametrized and non-sparse, have a higher computational cost. This motivates the search for optimal
sparse realizations [1]. Interestingly, the additional compu-
tational cost of the δ-realizations is small compared to the
very large decrease in sensitivities. This supports the view
that δ-realizations are generally superior [4].

VI. CONCLUSIONS

The Specialized Implicit Form provides a general frame-
work for the analysis and design of digital controller im-
plementations with minimal finite wordlength effects. This
paper has presented the closed-loop sensitivity measure in
this context. In this sense, it generalizes the results of [1].
Moreover, the present development applies for both fixed-
point and floating-point arithmetic, enabling one to analyze
more precisely the effects of quantization and rounding on
the parameters of a digital controller implementation.

As shown on two examples, the measure is computation-
ally tractable and hence amenable to solving the problem
of finding a good realization with regard to the parametric
sensitivity FWL effect. Furthermore, the sensitivity to the
rounding of each individual parameter can be easily obtained,
which may be very useful from a methodological point of
view.

Although the computational cost has been given for each
presented realization, they have not been discussed. The way
to manage the compromise between the sensitivity measure,
the computational effort, the rounding noise and the risk of
overflow will be considered in a future paper.

REFERENCES

[1] M. Gevers and G. Li, *Parametrizations in Control, Estimation and Filtering Problems*. Springer-Verlag, 1993.

[2] R. Istepanian and J. Whidborne, Eds., *Digital Controller Implementation and Fragility: A Modern Perspective*. London, UK: Springer-Verlag, Sept. 2001.

[3] J. Whidborne, R. Istepanian, and J. Wu, "Reduction of controller fragility by pole sensitivity minimization," *IEEE Trans. Automatic Control*, vol. 46, pp. 320–325, 2001.

[4] R. Goodall, "Perspectives on processing for real-time control," *Annual Reviews in Control*, vol. 25, pp. 123–131, 2001.

[5] L. Thiele, "Design of sensitivity and round-off noise optimal state-space discrete systems," *Int. J. Circuit Theory Appl.*, vol. 12, pp. 39–46, 1984.

[6] T. Hilaire, P. Chevrel, and Y. Trinquet, "Implicit state-space representation : a unifying framework for FWL implementation of LTI systems," in *Proc. of the 16th IFAC World Congress*, P. Piztek, Ed. Elsevier, July 2005.

[7] ——, "Designing low parametric sensitivity FWL realizations of LTI controllers/filters within the implicit state-space framework," in *Proc. of the 44th IEEE Conference on Decision and Control and the European Control Conference (CDC-ECC'05)*, December 2005, pp. 5192–5197.

[8] T. Hilaire, "Analyse et synthèse de l'implémentation de lois de contrôle-commande en précision finie (étude dans le cadre des applications automobiles sur calculateur embarquée)," Ph.D. dissertation, Université de Nantes, June 2006.

[9] J. Aplevich, *Implicit Linear Systems*. Springer-Verlag, 1991.

[10] T. Hilaire, P. Chevrel, and J. Whidborne, "A unifying framework for finite wordlength realizations," *IEEE Trans. on Circuits and Systems (accepted)*, 2007.

[11] J. Wu, S. Chen, and J. Chu, "Comparative study on finite-precision controller realizations in different representation schemes," *9th Annual Conf. Chinese Automation and Computing Society in UK (Luton, UK)*, September 2003.

[12] I. Fialho and T. Georgiou, "On stability and performance of sampled-data systems subject to wordlength constraint'," in *IEEE Trans. Automatic Control*, vol. 39, no. 12, December 1994, pp. 2476–2481.

[13] J. Whidborne, J. Wu, and R. Istepanian, "Finite word length stability issues in an $\ell_1$ framework," *Int. J. Control*, vol. 73, no. 2, pp. 166–176, 2000.

[14] G. Li, "On the structure of digital controllers with finite word length consideration," in *IEEE Trans. on Autom. Control*, vol. 43, no. 5, May 1998, pp. 689–693.

[15] J. Wu, S. Chen, G. Li, R. Istepanian, and J. Chu, "An improved closed-loop stability related measure for finite-precision digital controller realizations," *IEEE Trans. Automatic Control*, vol. 46, no. 7, pp. 1162–1166, 2001.

[16] W.-S. Yu and H.-J. Ko, "Improved eigenvalue sensitivity for finite-precision digital controller realisations via orthogonal Hermitian transform," *IEE Proc. Control Theory and Appl.*, vol. 150, no. 4, pp. 365–375, 2003.

[17] K. Zhou, J. Doyle, and K. Glover, *Robust and Optimal Control*. Prentice-Hall, 1996.

[18] T. Hilaire, P. Chevrel, and J.-P. Clauzel, "Low parametric sensitivity realization design for fwl implementation of mimo controllers : Theory and application to the active control of vehicle longitudinal oscillations," in *Proc. of Control Applications of Optimisation CAO'06*, April 2006.

[19] M. Ikeda, D. Šiljak, and D. White, "An inclusion principle for dynamic systems," *IEEE Trans. Automatic Control*, vol. 29, no. 3, pp. 244–249, March 1984.

[20] D. Šiljak, *Decentralized Control of Complex Systems*. Academic Press, 1991.

[21] L. Ingber, "Adaptive simulated annealing (ASA): Lessons learned," *Control and Cybernetics*, vol. 25, no. 1, pp. 33–54, 1996.

[22] S. Chen and B. Luk, "Adaptive simulated annealing for optimization in signal processing applications," *Signal Processing*, vol. 79, pp. 117–128, 1999.

[23] I. Njabeleke, R. Pannett, P. Chawdhry, and C. Burrows, "Design of h-infiny loop-shaping controllers for fluid power systems," *IEE Colloquium Robust Control - Theory, Software and Applications*, 1997.

[24] D. Lefebvre, P. Chevrel, and S. Richard, "An h-infinity based control design methodology dedicated to the active control of longitudinal oscillations," *IEEE Trans. on Control Systems Technology*, vol. 11, no. 6, pp. 948–956, November 2003.

[25] D. Alazard and P. Apkarian, "Exact observer-based structures for arbitrary compensators," *International Journal of Robust and Non-Linear Control*, no. 9, pp. 101–118, 1999.
APPENDIX

$$A_p = \begin{pmatrix} 9.9988e-1 & 1.9432e-5 & 5.9320e-5 & -6.2286e-5 \\ -4.9631e-7 & 2.3577e-2 & 2.3709e-5 & 2.3672e-5 \\ -1.5151e-3 & 2.3709e-2 & 2.3751e-5 & 2.3898e-5 \\ 1.5908e-3 & 2.3672e-2 & 2.3898e-5 & 2.3667e-5 \end{pmatrix} B_p = \begin{pmatrix} 3.0504e-3 & -1.2373e-2 & -1.2375e-2 & -8.8703e-2 \end{pmatrix} C_p = \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \end{pmatrix}^\top \quad (48)$$

$$Z_0 = \begin{pmatrix} 3.325e+0 & -3.982e+0 & 1.987e+0 & -3.307e-1 & 1 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ -1.611e-3 & 3.758e-3 & -2.684e-3 & 5.366e-4 & -8.084e-4 \end{pmatrix}, \quad \frac{\delta \hat{H}}{\delta Z} \Big|_{Z_0} = \begin{pmatrix} 6.291e+7 & 6.291e+7 & 6.291e+7 & 6.291e+7 & 6.065e-2 \\ 1.464e+8 & 1.464e+8 & 1.464e+8 & 1.464e+8 & 1.413e-1 \\ 1.044e+8 & 1.044e+8 & 1.044e+8 & 1.044e+8 & 1.008e-1 \\ 2.085e+7 & 2.085e+7 & 2.085e+7 & 2.085e+7 & 2.016e-2 \\ 1.536e+8 & 1.536e+8 & 1.536e+8 & 1.536e+8 & 2.315e-1 \end{pmatrix} \quad (49)$$

$$Z_1 = \begin{pmatrix} 7.467e-1 & -1.440e-1 & -1.573e-1 & -1.408e-2 & -1.954e-4 \\ -5.086e-1 & 6.942e-1 & -3.288e-1 & -2.688e-2 & 3.935e-3 \\ -1.657e-1 & -9.209e-2 & 8.991e-1 & -9.481e-3 & -3.448e-3 \\ -2.470e-1 & -1.388e-1 & -1.501e-1 & 9.855e-1 & 1.496e-3 \\ 1.430e+0 & -7.149e-1 & -3.335e-1 & 2.217e-1 & -8.084e-4 \end{pmatrix}, \quad \frac{\delta \hat{H}}{\delta Z} \Big|_{Z_1} = \begin{pmatrix} 1.353e+2 & 2.130e+2 & 3.249e+2 & 8.355e+2 & 4.577e+1 \\ 6.131e+1 & 9.905e+1 & 1.492e+2 & 3.784e+2 & 2.194e+1 \\ 7.568e+1 & 9.976e+1 & 1.672e+2 & 4.685e+2 & 1.473e+1 \\ 2.167e+1 & 3.359e+1 & 5.162e+1 & 1.338e+2 & 7.081e+0 \\ 6.393e-1 & 1.037e+0 & 1.559e+0 & 3.946e+0 & 2.315e-1 \end{pmatrix} \quad (50)$$

$$Z_3 = \begin{pmatrix} -1 & 0 & 0 & 0 & -8.316e+0 & 1.114e+1 & 3.108e+0 & -4.865e+0 & 3.683e-2 \\ 0 & -1 & 0 & 0 & 6.036e+0 & -8.458e+0 & -2.466e+0 & 3.662e+0 & 2.341e-2 \\ 0 & 0 & -1 & 0 & -8.472e+0 & 1.157e+1 & 3.284e+0 & -5.042e+0 & 4.442e-3 \\ 0 & 0 & 0 & -1 & -1.347e+1 & 1.863e+1 & 5.358e+0 & -8.094e+0 & -4.777e-3 \\ 3.125e-2 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 3.125e-2 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 3.125e-2 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 3.125e-2 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & -1.248e+1 & 1.525e+1 & 4.116e+0 & -6.892e+0 & -8.084e-4 \end{pmatrix} \quad (51)$$

$$\frac{\delta \hat{H}}{\delta Z} \Big|_{Z_3} = \begin{pmatrix} 5.131e-3 & 8.472e-3 & 7.613e-3 & 1.031e-2 & 2.633e-1 & 6.329e-1 & 5.877e-1 & 8.116e-1 & 1.423e+0 \\ 3.931e-3 & 7.367e-3 & 6.688e-3 & 9.210e-3 & 3.071e-1 & 1.103e+0 & 1.034e+0 & 1.457e+0 & 1.260e+0 \\ 8.230e-3 & 1.553e-2 & 1.418e-2 & 1.934e-2 & 5.398e-1 & 1.579e+0 & 1.475e+0 & 2.060e+0 & 2.659e+0 \\ 4.606e-3 & 9.847e-3 & 9.082e-3 & 1.247e-2 & 3.818e-1 & 1.265e+0 & 1.185e+0 & 1.664e+0 & 1.708e+0 \\ 1.642e-1 & 2.711e-1 & 2.436e-1 & 3.300e-1 & 8.426e+0 & 2.025e+1 & 1.881e+1 & 2.597e+1 & 4.554e+1 \\ 1.258e-1 & 2.357e-1 & 2.140e-1 & 2.947e-1 & 9.826e+0 & 3.531e+1 & 3.309e+1 & 4.664e+1 & 4.033e+1 \\ 2.634e-1 & 4.970e-1 & 4.538e-1 & 6.189e-1 & 1.728e+1 & 5.053e+1 & 4.720e+1 & 6.593e+1 & 8.508e+1 \\ 1.474e-1 & 3.151e-1 & 2.906e-1 & 3.990e-1 & 1.222e+1 & 4.048e+1 & 3.790e+1 & 5.324e+1 & 5.466e+1 \\ 1.079e-3 & 1.442e-3 & 1.243e-3 & 1.669e-3 & 4.075e-2 & 8.483e-2 & 7.809e-2 & 1.071e-1 & 2.315e-1 \end{pmatrix} \quad (52)$$

$$A_p = \begin{pmatrix} 8.384e-1 & 1.600e-1 & -3.294e-1 & -4.833e-2 & 0 & 0 & 0 & 0 & 0 & 0 \\ -3.927e-1 & 7.144e-1 & 5.040e-2 & -8.245e-3 & 0 & 0 & 0 & 0 & 0 & 0 \\ -1.566e-1 & -6.105e-1 & 3.683e-2 & 4.195e-1 & 0 & 0 & 0 & 0 & 0 & 0 \\ -1.444e-1 & 1.772e-1 & -6.798e-1 & 6.508e-1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1.525e-1 & 1.512e-1 & 4.030e-1 & 3.898e-1 & 9.773e-1 & 1.037e-2 & -6.170e-2 & 0 & 0 & 0 \\ 2.768e-4 & 2.170e-4 & 5.783e-4 & 5.594e-4 & 2.837e-3 & 9.971e-1 & 1.698e-2 & 0 & 0 & 0 \\ 3.238e-2 & 2.539e-2 & 6.767e-2 & 6.545e-2 & 3.320e-1 & -3.341e-1 & 9.868e-1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1.000e+0 & -1.000e-10 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1.000e-2 & 1.000e+0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 9.417e-1 \end{pmatrix} \quad (53)$$

$$B_p = \begin{pmatrix} -4.007e+0 & 0 & 0 & 0 & -2.372e+0 & 0 & 0 & 0 & 0 & 0 \\ -5.769e+0 & 0 & 0 & 0 & -2.540e+0 & 0 & 0 & 0 & 0 & 0 \\ -6.522e+0 & 0 & 0 & 0 & -1.210e-1 & 0 & 0 & 0 & 0 & 0 \\ 2.490e+0 & 0 & 0 & 0 & -1.565e-4 & 0 & 0 & 0 & 0 & 0 \\ 8.562e-1 & 0 & 0 & 0 & -6.245e-2 & 0 & 0 & 0 & 0 & 0 \\ 1.229e-3 & 0 & 0 & 0 & 1.151e+0 & 0 & 0 & 0 & 0 & 0 \\ 1.438e-1 & 0 & 0 & 0 & 4.083e-2 & 0 & 0 & 0 & 0 & 0 \\ 1.000e+0 & 0 & 0 & 0 & 2.255e-1 & 0 & 0 & 0 & 0 & 0 \\ 5.000e-3 & 0 & 0 & 0 & -1.528e-2 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -9.720e-4 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, C_p = \begin{pmatrix} 9.209e-3 & 0 & 0 & 0 & -2.372e+0 & 0 & 0 & 0 & 0 & 0 \\ 7.221e-3 & 0 & 0 & 0 & -2.540e+0 & 0 & 0 & 0 & 0 & 0 \\ 1.924e-2 & 0 & 0 & 0 & -1.210e-1 & 0 & 0 & 0 & 0 & 0 \\ 1.861e-2 & 0 & 0 & 0 & -1.565e-4 & 0 & 0 & 0 & 0 & 0 \\ 9.441e-2 & 0 & 0 & 0 & -6.245e-2 & 0 & 0 & 0 & 0 & 0 \\ 4.953e-4 & 0 & 0 & 0 & 1.151e+0 & 0 & 0 & 0 & 0 & 0 \\ -2.946e-3 & 0 & 0 & 0 & 4.083e-2 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 2.255e-1 & 0 & 0 & 0 & 0 & 0 \\ -3.495e-1 & 0 & 0 & 0 & -1.528e-2 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, C = \begin{pmatrix} -2.372e-2 & 0 & 0 & 0 & 2.540e-2 & 0 & 0 & 0 & 0 & 0 \\ -1.210e-3 & 0 & 0 & 0 & 1.210e-3 & 0 & 0 & 0 & 0 & 0 \\ -1.565e-6 & 0 & 0 & 0 & 6.245e-4 & 0 & 0 & 0 & 0 & 0 \\ 1.151e-2 & 0 & 0 & 0 & 4.083e-4 & 0 & 0 & 0 & 0 & 0 \\ -2.255e-3 & 0 & 0 & 0 & 1.528e-4 & 0 & 0 & 0 & 0 & 0 \\ 9.720e-6 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, D = -2.140e-1 \quad (54)$$

$$A = \begin{pmatrix} 8.195e-1 & 2.812e-1 & -3.317e-2 & 2.699e-2 & -1.649e-1 & 1.318e-1 & 1.059e-2 & -6.733e-2 & 1.750e-3 & 6.525e-5 \\ -2.812e-1 & -4.817e-1 & -1.668e-1 & 8.654e-2 & -5.403e-1 & 1.469e-1 & 1.837e-2 & -1.211e-1 & 1.942e-3 & 2.134e-5 \\ 3.317e-2 & -1.668e-1 & 9.749e-1 & 1.696e-2 & -9.104e-2 & 7.638e-2 & 3.357e-3 & -2.006e-2 & 8.441e-4 & 4.548e-5 \\ 2.699e-2 & -8.654e-2 & -1.696e-2 & 9.601e-1 & 2.528e-1 & 5.956e-2 & 1.654e-3 & -9.085e-3 & 6.046e-4 & 3.843e-5 \\ 1.649e-1 & -5.403e-1 & -9.104e-2 & -2.528e-1 & 6.022e-1 & 3.888e-1 & 1.150e-2 & -6.420e-2 & 3.945e-3 & 2.454e-4 \\ 1.318e-1 & -1.469e-1 & -7.638e-2 & 5.956e-2 & -3.888e-1 & 4.664e-1 & -6.206e-2 & 4.224e-1 & -8.490e-4 & 3.703e-4 \\ 1.059e-2 & -1.837e-2 & -3.357e-3 & 1.654e-3 & -1.150e-2 & -6.206e-2 & 9.832e-1 & 1.258e-1 & 7.737e-3 & 6.392e-4 \\ 6.733e-2 & -1.211e-1 & -2.006e-2 & 9.085e-3 & -6.420e-2 & -4.224e-1 & -1.258e-1 & -4.483e-2 & 7.258e-2 & 5.631e-3 \\ -1.750e-3 & 1.942e-3 & 8.441e-4 & -6.046e-4 & 3.945e-3 & 8.490e-4 & 7.258e-2 & 9.838e-1 & -2.474e-3 & 9.418e-1 \end{pmatrix} \quad (55)$$

$$Z = \begin{pmatrix} I & 0 & -C_p & I \\ -Q & I & -K_c & 0 \\ -K_f & -B_p & A_p & 0 \\ 0 & -I & 0 & 0 \end{pmatrix} \quad (56)$$