After completing this lesson, students should be able to:

• represent a multivariable feedback/feedforward controller by one matrix parameter and one regressor vector;
• convert matrix equations into standard linear regressions by vectorization and Kronecker products;
• determine whether a desired reference model is exactly matchable by a chosen controller structure;
• impose known sparsity, decentralization, symmetry, or basis constraints without introducing redundant adaptive parameters;
• distinguish a controller that is linear in adjustable parameters from a controller whose internal dynamics make the parameter dependence nonlinear;
• test whether the chosen regressor can distinguish the adjustable multivariable parameters.

Chapter 23, Lesson 1 introduced parametric uncertainty in MIMO plants. Here, the problem is not yet to prove stability of a particular adaptive law. The immediate problem is more basic: what quantities should be adapted, and how should they be arranged? A poor parameterization may violate known interconnection structure, create unidentifiable directions, or make the controller nonlinear in its adjustable parameters.

Consider a continuous-time plant with \( n \) states, \( m \) control inputs, and \( p \) measured outputs:

\[ \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t), \qquad \mathbf{y}(t)=\mathbf{C}\mathbf{x}(t), \]

where \( \mathbf{x}\in\mathbb{R}^{n} \), \( \mathbf{u}\in\mathbb{R}^{m} \), and \( \mathbf{y}\in\mathbb{R}^{p} \). For state tracking, choose a stable reference model

\[ \dot{\mathbf{x} }_m(t) =\mathbf{A}_m\mathbf{x}_m(t)+\mathbf{B}_m\mathbf{r}(t), \qquad \mathbf{y}_m(t)=\mathbf{C}_m\mathbf{x}_m(t), \]

with command \( \mathbf{r}\in\mathbb{R}^{q} \). The most direct static full-state controller parameterization is

\[ \mathbf{u}(t) =\mathbf{K}_x(t)\mathbf{x}(t) +\mathbf{K}_r(t)\mathbf{r}(t), \]

where

\[ \mathbf{K}_x(t)\in\mathbb{R}^{m\times n}, \qquad \mathbf{K}_r(t)\in\mathbb{R}^{m\times q}. \]

Stack the two gain matrices horizontally and the signals vertically:

\[ \boldsymbol{\Theta}(t) = \begin{bmatrix} \mathbf{K}_x(t) & \mathbf{K}_r(t) \end{bmatrix} \in\mathbb{R}^{m\times(n+q)}, \qquad \boldsymbol{\phi}(t) = \begin{bmatrix} \mathbf{x}(t)\\ \mathbf{r}(t) \end{bmatrix} \in\mathbb{R}^{n+q}. \]

The controller then becomes the compact matrix regression

\[ \boxed{\mathbf{u}(t)=\boldsymbol{\Theta}(t)\boldsymbol{\phi}(t)}. \]

A dense parameterization contains \( m(n+q) \) independently adjustable scalar parameters. Every input channel may depend on every state and every reference component.

flowchart TD
  R["Command vector r"] --> RM["Fixed reference model"]
  RM --> XM["Desired state xm"]
  R --> REG["Regressor phi = stack(x, r)"]
  X["Measured plant state x"] --> REG
  REG --> CTRL["Controller u = Theta phi"]
  CTRL --> P["MIMO plant"]
  P --> X
  X --> ERR["Tracking error e = x - xm"]
  XM --> ERR
  ERR --> ADAPT["Matrix parameter update"]
  REG --> ADAPT
  ADAPT --> CTRL
        

Most adaptive-law derivations are written for a parameter vector. Define the column-stacked parameter vector

\[ \boldsymbol{\theta}(t) =\operatorname{vec}\!\left(\boldsymbol{\Theta}(t)\right) \in\mathbb{R}^{m(n+q)}. \]

The matrix controller can be written as a conventional vector regression:

\[ \boxed{ \mathbf{u}(t) = \left( \boldsymbol{\phi}^{T}(t)\otimes\mathbf{I}_{m} \right)\boldsymbol{\theta}(t) } = \boldsymbol{\Psi}(t)\boldsymbol{\theta}(t), \]

where

\[ \boldsymbol{\Psi}(t) = \boldsymbol{\phi}^{T}(t)\otimes\mathbf{I}_{m} \in\mathbb{R}^{m\times m(n+q)}. \]

3.1 Proof of the identity

Let \( \boldsymbol{\Theta} =[\boldsymbol{\theta}_1\ \cdots\ \boldsymbol{\theta}_{n+q}] \), where each column \( \boldsymbol{\theta}_j\in\mathbb{R}^{m} \). Then

\[ \boldsymbol{\Theta}\boldsymbol{\phi} = \sum_{j=1}^{n+q} \boldsymbol{\theta}_{j}\phi_j. \]

Column-wise vectorization gives

\[ \operatorname{vec}(\boldsymbol{\Theta}) = \begin{bmatrix} \boldsymbol{\theta}_{1}\\ \vdots\\ \boldsymbol{\theta}_{n+q} \end{bmatrix}, \qquad \boldsymbol{\phi}^{T}\otimes\mathbf{I}_{m} = \begin{bmatrix} \phi_1\mathbf{I}_{m} & \cdots & \phi_{n+q}\mathbf{I}_{m} \end{bmatrix}. \]

Multiplication therefore yields

\[ \left( \boldsymbol{\phi}^{T}\otimes\mathbf{I}_{m} \right) \operatorname{vec}(\boldsymbol{\Theta}) = \sum_{j=1}^{n+q}\phi_j\boldsymbol{\theta}_j = \boldsymbol{\Theta}\boldsymbol{\phi}. \]

This identity is central: it permits matrix gains to be adapted either directly in matrix form or through a single scalar-parameter vector. The two descriptions are algebraically equivalent as long as the same vectorization convention is used.

3.2 Matrix gradient and vector gradient

For the instantaneous quadratic control-output loss

\[ J(\boldsymbol{\Theta}) = \frac{1}{2} \left\| \mathbf{u}_d-\boldsymbol{\Theta}\boldsymbol{\phi} \right\|_2^2, \]

define

\[ \boldsymbol{\varepsilon}_u = \boldsymbol{\Theta}\boldsymbol{\phi}-\mathbf{u}_d. \]

The Frobenius gradient is

\[ \nabla_{\boldsymbol{\Theta} }J = \boldsymbol{\varepsilon}_u\boldsymbol{\phi}^{T}. \]

Vectorizing this relation gives the usual regression gradient

\[ \nabla_{\boldsymbol{\theta} }J = \boldsymbol{\Psi}^{T}\boldsymbol{\varepsilon}_u. \]

Thus, an outer product in matrix coordinates is the same operation as a regressor-transpose multiplication in vector coordinates.

Suppose the plant state is available and the controller has ideal constant gains \( \mathbf{K}_x^{*} \) and \( \mathbf{K}_r^{*} \). Substitution into the plant gives

\[ \dot{\mathbf{x} } = \left( \mathbf{A}+\mathbf{B}\mathbf{K}_x^{*} \right)\mathbf{x} + \mathbf{B}\mathbf{K}_r^{*}\mathbf{r}. \]

Exact state-model matching requires the closed-loop plant dynamics to equal the reference-model dynamics:

\[ \boxed{ \mathbf{A}+\mathbf{B}\mathbf{K}_x^{*}=\mathbf{A}_m, \qquad \mathbf{B}\mathbf{K}_r^{*}=\mathbf{B}_m. } \]

These are algebraic feasibility conditions. An adaptive law cannot overcome failure of the selected controller structure to satisfy them.

4.1 Existence theorem

The equation \( \mathbf{B}\mathbf{K}=\mathbf{D} \) is solvable if and only if every column of \( \mathbf{D} \) lies in the column space of \( \mathbf{B} \). Equivalently,

\[ \left( \mathbf{I}_{n}-\mathbf{B}\mathbf{B}^{\dagger} \right)\mathbf{D} = \mathbf{0}, \]

where \( \mathbf{B}^{\dagger} \) is the Moore-Penrose pseudoinverse.

Proof. If a solution exists, then

\[ \mathbf{D}=\mathbf{B}\mathbf{K}, \]

so every column of \( \mathbf{D} \) is a linear combination of columns of \( \mathbf{B} \). Conversely, if \( \mathbf{D} \) belongs to the column space of \( \mathbf{B} \), then

\[ \mathbf{B}\mathbf{B}^{\dagger}\mathbf{D} = \mathbf{D}, \]

and \( \mathbf{K}=\mathbf{B}^{\dagger}\mathbf{D} \) is one solution.

Applying the theorem to the matching equations gives the feasibility tests

\[ \left( \mathbf{I}_{n}-\mathbf{B}\mathbf{B}^{\dagger} \right) \left( \mathbf{A}_m-\mathbf{A} \right) = \mathbf{0}, \]

\[ \left( \mathbf{I}_{n}-\mathbf{B}\mathbf{B}^{\dagger} \right) \mathbf{B}_m = \mathbf{0}. \]

4.2 Complete family of ideal gains

If the matching equations are feasible, all solutions are

\[ \mathbf{K}_x^{*} = \mathbf{B}^{\dagger} \left( \mathbf{A}_m-\mathbf{A} \right) + \left( \mathbf{I}_{m}-\mathbf{B}^{\dagger}\mathbf{B} \right)\mathbf{Z}_x, \]

\[ \mathbf{K}_r^{*} = \mathbf{B}^{\dagger}\mathbf{B}_m + \left( \mathbf{I}_{m}-\mathbf{B}^{\dagger}\mathbf{B} \right)\mathbf{Z}_r, \]

where \( \mathbf{Z}_x \) and \( \mathbf{Z}_r \) are arbitrary matrices of compatible dimensions. The null-space terms do not change \( \mathbf{B}\mathbf{K}_x^{*} \) or \( \mathbf{B}\mathbf{K}_r^{*} \).

If \( \mathbf{B}\in\mathbb{R}^{n\times n} \) is nonsingular, the ideal gains are unique:

\[ \mathbf{K}_x^{*} = \mathbf{B}^{-1} \left( \mathbf{A}_m-\mathbf{A} \right), \qquad \mathbf{K}_r^{*} = \mathbf{B}^{-1}\mathbf{B}_m. \]

For an overactuated plant, the matching gains are generally nonunique. For an underactuated plant, exact matching is generally restrictive because the columns of the desired dynamic correction must lie in the input subspace.

The dense matrix \( \boldsymbol{\Theta} \) is maximally flexible, but it may ignore known architecture. A structured parameterization should preserve genuine prior information while retaining enough freedom to satisfy the matching equations.

5.1 Affine basis parameterization

Let \( \mathbf{F}_1,\ldots,\mathbf{F}_{s} \) be fixed controller basis matrices. Write

\[ \boldsymbol{\Theta}(\boldsymbol{\alpha}) = \boldsymbol{\Theta}_0 + \sum_{i=1}^{s}\alpha_i\mathbf{F}_i, \qquad \boldsymbol{\alpha}\in\mathbb{R}^{s}. \]

Then the control input is linear in \( \boldsymbol{\alpha} \):

\[ \mathbf{u} = \boldsymbol{\Theta}_0\boldsymbol{\phi} + \begin{bmatrix} \mathbf{F}_1\boldsymbol{\phi} & \cdots & \mathbf{F}_{s}\boldsymbol{\phi} \end{bmatrix} \boldsymbol{\alpha}. \]

Define the reduced regressor

\[ \boldsymbol{\Psi}_{\alpha}(\boldsymbol{\phi}) = \begin{bmatrix} \mathbf{F}_1\boldsymbol{\phi} & \cdots & \mathbf{F}_{s}\boldsymbol{\phi} \end{bmatrix} \in\mathbb{R}^{m\times s}. \]

The number of adapted scalars falls from \( m(n+q) \) to \( s \).

5.2 Selection-matrix form

Vectorize the affine basis:

\[ \operatorname{vec} \left( \boldsymbol{\Theta}(\boldsymbol{\alpha}) \right) = \boldsymbol{\theta}_0+\mathbf{S}\boldsymbol{\alpha}, \]

where

\[ \mathbf{S} = \begin{bmatrix} \operatorname{vec}(\mathbf{F}_1) & \cdots & \operatorname{vec}(\mathbf{F}_{s}) \end{bmatrix}. \]

The corresponding reduced vector regression is

\[ \mathbf{u} = \boldsymbol{\Psi}\boldsymbol{\theta}_0 + \boldsymbol{\Psi}\mathbf{S}\boldsymbol{\alpha}. \]

5.3 Masked sparsity and decentralized structures

Let \( \mathbf{M}\in\{0,1\}^{m\times(n+q)} \) be a binary mask. A masked controller is

\[ \boldsymbol{\Theta} = \boldsymbol{\Theta}_0 + \mathbf{M}\odot\mathbf{W}, \]

where \( \odot \) is element-wise multiplication. A zero in \( \mathbf{M} \) fixes the corresponding controller entry. This is appropriate when communication or sensing constraints genuinely prohibit a channel.

For two weakly coupled subsystems, a block-diagonal feedback matrix has the form

\[ \mathbf{K}_x = \begin{bmatrix} \mathbf{K}_{x,1} & \mathbf{0}\\ \mathbf{0} & \mathbf{K}_{x,2} \end{bmatrix}. \]

Such a decentralized parameterization reduces adaptation complexity, but it can make exact matching impossible when important plant interactions require off-diagonal compensation.

5.4 Shared-parameter and symmetry constraints

Suppose physical symmetry implies repeated gains. For example,

\[ \mathbf{K}_x(\alpha_1,\alpha_2) = \begin{bmatrix} \alpha_1 & \alpha_2\\ \alpha_2 & \alpha_1 \end{bmatrix}. \]

This structure uses two scalars instead of four. Its basis matrices are

\[ \mathbf{F}_1 = \begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix}, \qquad \mathbf{F}_2 = \begin{bmatrix} 0 & 1\\ 1 & 0 \end{bmatrix}. \]

5.5 Why low-rank factorization is different

A low-rank controller can be written as

\[ \boldsymbol{\Theta} = \mathbf{U}\mathbf{V}^{T}. \]

Although this may reduce storage, the control input is bilinear in the adjustable factors \( \mathbf{U} \) and \( \mathbf{V} \). Therefore, standard linear parameter-adaptation arguments do not apply directly. If linear parameterization is required, a fixed basis with adaptive scalar coefficients is usually preferable.

flowchart TD
  A["Start from matching and implementation requirements"] --> B["Is every cross-channel gain physically available?"]
  B -->|yes| C["Use full dense matrix as baseline"]
  B -->|no| D["Impose a sparsity or block mask"]
  C --> E["Is reliable symmetry or repeated-gain information known?"]
  D --> E
  E -->|yes| F["Build an affine basis with shared parameters"]
  E -->|no| G["Keep independent allowed entries"]
  F --> H["Check matching equations inside the chosen subspace"]
  G --> H
  H -->|feasible| I["Form reduced Kronecker regressor"]
  H -->|not feasible| J["Enlarge structure or redesign reference model"]
        

Output-feedback controllers often require dynamic compensation. A safe way to preserve linearity in adjustable parameters is to generate regressors through fixed stable filters and adapt only the output weights.

Let \( \boldsymbol{\omega}_y \) and \( \boldsymbol{\omega}_r \) satisfy

\[ \dot{\boldsymbol{\omega} }_y = \mathbf{F}_y\boldsymbol{\omega}_y + \mathbf{G}_y\mathbf{y}, \qquad \dot{\boldsymbol{\omega} }_r = \mathbf{F}_r\boldsymbol{\omega}_r + \mathbf{G}_r\mathbf{r}, \]

where \( \mathbf{F}_y \) and \( \mathbf{F}_r \) are fixed Hurwitz matrices. Form

\[ \boldsymbol{\omega} = \begin{bmatrix} \boldsymbol{\omega}_y\\ \boldsymbol{\omega}_r\\ \mathbf{y}\\ \mathbf{r} \end{bmatrix}, \qquad \mathbf{u} = \boldsymbol{\Theta}_{\omega}\boldsymbol{\omega}. \]

Because the filter dynamics are fixed, \( \boldsymbol{\omega}(t) \) is a measurable regressor independent of the current adjustable matrix \( \boldsymbol{\Theta}_{\omega} \). The control law is therefore linear in its adjustable parameters.

6.1 A parameter-dependent controller state is generally nonlinear

Consider instead

\[ \dot{\boldsymbol{\xi} } = \mathbf{A}_c(\boldsymbol{\eta})\boldsymbol{\xi} + \mathbf{B}_c(\boldsymbol{\eta})\mathbf{y}, \qquad \mathbf{u} = \mathbf{C}_c(\boldsymbol{\eta})\boldsymbol{\xi}. \]

Even if each controller matrix is affine in \( \boldsymbol{\eta} \), the controller state \( \boldsymbol{\xi}(t) \) also depends on the entire history of \( \boldsymbol{\eta} \). Consequently, \( \mathbf{u}(t) \) is not generally a linear regression in the current parameter vector. This distinction is essential when importing gradient or Lyapunov adaptive laws.

7.1 Keep the learning target fixed

The reference model defines the desired closed-loop behavior. In the basic MRAC architecture, its matrices are selected in advance and are not adapted:

\[ \dot{\mathbf{x} }_m = \mathbf{A}_m\mathbf{x}_m + \mathbf{B}_m\mathbf{r}. \]

Adapting the controller while simultaneously moving the reference model changes the target and complicates interpretation of the tracking error. Reference-model modification is a separate advanced design choice, not a substitute for a feasible controller parameterization.

7.2 Coordinate transformation of the gain matrix

Let \( \mathbf{x}=\mathbf{T}\mathbf{z} \), where \( \mathbf{T} \) is nonsingular. Then

\[ \mathbf{u} = \mathbf{K}_x\mathbf{x} = \mathbf{K}_x\mathbf{T}\mathbf{z} = \mathbf{K}_z\mathbf{z}, \qquad \mathbf{K}_z=\mathbf{K}_x\mathbf{T}. \]

A diagonal \( \mathbf{K}_x \) does not generally remain diagonal after an arbitrary state transformation. Therefore, sparsity and decentralization constraints are coordinate-dependent. They should be imposed in coordinates with a physical interpretation, such as local subsystem states or actual communication channels.

7.3 Reference command scaling

If a nonsingular command transformation \( \mathbf{r}=\mathbf{T}_r\boldsymbol{\rho} \) is used, then

\[ \mathbf{K}_r\mathbf{r} = \mathbf{K}_r\mathbf{T}_r\boldsymbol{\rho}, \qquad \mathbf{K}_{\rho} = \mathbf{K}_r\mathbf{T}_r. \]

The feedforward parameter matrix must be transformed consistently; otherwise, numerical parameter values may be compared across incompatible command coordinates.

Let \( \mathbf{s}(t)\in\mathbb{R}^{m} \) be an adaptation-driving error signal obtained from a tracking error, filtered error, or prediction error. A generic matrix update is

\[ \dot{\boldsymbol{\Theta} } = -\mathbf{\Gamma}_u \mathbf{s}\boldsymbol{\phi}^{T} \mathbf{\Gamma}_{\phi}, \]

where \( \mathbf{\Gamma}_u\in\mathbb{R}^{m\times m} \) and \( \mathbf{\Gamma}_{\phi}\in \mathbb{R}^{(n+q)\times(n+q)} \) are positive-definite adaptation-gain matrices. The left gain scales input-channel update directions, and the right gain scales regressor directions.

In vector form,

\[ \dot{\boldsymbol{\theta} } = - \left( \mathbf{\Gamma}_{\phi}^{T}\otimes\mathbf{\Gamma}_u \right) \boldsymbol{\Psi}^{T}\mathbf{s}. \]

8.1 Derivation by vectorization

Using

\[ \operatorname{vec} \left( \mathbf{A}\mathbf{X}\mathbf{B} \right) = \left( \mathbf{B}^{T}\otimes\mathbf{A} \right) \operatorname{vec}(\mathbf{X}), \]

we obtain

\[ \operatorname{vec} \left( \mathbf{\Gamma}_u \mathbf{s}\boldsymbol{\phi}^{T} \mathbf{\Gamma}_{\phi} \right) = \left( \mathbf{\Gamma}_{\phi}^{T}\otimes\mathbf{\Gamma}_u \right) \operatorname{vec} \left( \mathbf{s}\boldsymbol{\phi}^{T} \right). \]

Furthermore,

\[ \operatorname{vec} \left( \mathbf{s}\boldsymbol{\phi}^{T} \right) = \left( \boldsymbol{\phi}\otimes\mathbf{I}_{m} \right)\mathbf{s} = \boldsymbol{\Psi}^{T}\mathbf{s}. \]

8.2 Normalization

To reduce sensitivity to large state or command magnitudes, use

\[ \nu(t) = 1+\boldsymbol{\phi}^{T}(t)\boldsymbol{\phi}(t), \qquad \dot{\boldsymbol{\Theta} } = -\frac{\gamma}{\nu(t)} \mathbf{s}\boldsymbol{\phi}^{T}. \]

The denominator does not change the update direction; it limits the instantaneous update magnitude when the regressor norm is large.

8.3 Projection and structural subspaces

If \( \boldsymbol{\Theta} =\boldsymbol{\Theta}_0+\operatorname{mat} (\mathbf{S}\boldsymbol{\alpha}) \), adapt only \( \boldsymbol{\alpha} \). A generic reduced update is

\[ \dot{\boldsymbol{\alpha} } = \operatorname{Proj}_{\mathcal{A} } \left[ \boldsymbol{\alpha}, -\mathbf{\Gamma}_{\alpha} \mathbf{S}^{T}\boldsymbol{\Psi}^{T}\mathbf{s} \right], \]

where \( \mathcal{A} \) is an admissible parameter set. This guarantees that adaptation cannot populate forbidden matrix entries.

The precise selection of \( \mathbf{s} \), signs, and gain matrices is tied to the MIMO error model and Lyapunov function. Those stability conditions are developed in Lesson 3.

A valid controller parameterization need not be identifiable from the available trajectories. Suppose

\[ \mathbf{u}(t) = \boldsymbol{\Psi}(t)\boldsymbol{\theta}. \]

Over an interval \( [t_0,t_1] \), define the information matrix

\[ \mathbf{G}_{\theta}(t_0,t_1) = \int_{t_0}^{t_1} \boldsymbol{\Psi}^{T}(\tau) \boldsymbol{\Psi}(\tau)\,d\tau. \]

Because

\[ \boldsymbol{\Psi}^{T}\boldsymbol{\Psi} = \left( \boldsymbol{\phi}\boldsymbol{\phi}^{T} \right) \otimes\mathbf{I}_{m}, \]

we have

\[ \mathbf{G}_{\theta} = \left( \int_{t_0}^{t_1} \boldsymbol{\phi}(\tau) \boldsymbol{\phi}^{T}(\tau)\,d\tau \right) \otimes\mathbf{I}_{m}. \]

9.1 Rank result

Using \( \operatorname{rank}(\mathbf{A}\otimes\mathbf{B}) =\operatorname{rank}(\mathbf{A}) \operatorname{rank}(\mathbf{B}) \),

\[ \operatorname{rank}(\mathbf{G}_{\theta}) = m\, \operatorname{rank} \left( \int_{t_0}^{t_1} \boldsymbol{\phi}\boldsymbol{\phi}^{T}\,d\tau \right). \]

Therefore, the dense matrix parameters can be distinguished only when the base regressor \( \boldsymbol{\phi} \) excites all \( n+q \) directions.

9.2 Structured information matrix

For \( \boldsymbol{\theta} =\boldsymbol{\theta}_0+\mathbf{S}\boldsymbol{\alpha} \), the reduced information matrix is

\[ \mathbf{G}_{\alpha} = \int_{t_0}^{t_1} \mathbf{S}^{T} \boldsymbol{\Psi}^{T}(\tau) \boldsymbol{\Psi}(\tau) \mathbf{S}\,d\tau. \]

If \( \mathbf{G}_{\alpha} \) is singular, distinct reduced parameter vectors generate the same control action along the observed trajectory. This may occur because:

• some states or commands remain zero;
• two regressor components are linearly dependent;
• the basis matrices themselves are linearly dependent;
• null-space controller directions do not affect the plant through \( \mathbf{B} \);
• a reference command is insufficiently rich for parameter convergence.

Tracking convergence and parameter convergence are different goals. The controller may track well while its individual matrix entries do not converge to a unique ideal value.

Consider

\[ \dot{\mathbf{x} } = \begin{bmatrix} -0.2 & 0.6\\ -0.4 & -0.3 \end{bmatrix} \mathbf{x} + \mathbf{u}, \]

and the desired reference model

\[ \dot{\mathbf{x} }_m = \begin{bmatrix} -1.5 & 0\\ 0 & -2 \end{bmatrix} \mathbf{x}_m + \mathbf{r}. \]

Here \( \mathbf{B}=\mathbf{I}_2 \) and \( \mathbf{B}_m=\mathbf{I}_2 \). The ideal gains are therefore

\[ \mathbf{K}_x^{*} = \mathbf{A}_m-\mathbf{A} = \begin{bmatrix} -1.3 & -0.6\\ 0.4 & -1.7 \end{bmatrix}, \qquad \mathbf{K}_r^{*} = \mathbf{I}_2. \]

With \( \boldsymbol{\phi} =[\mathbf{x}^{T}\ \mathbf{r}^{T}]^{T} \),

\[ \boldsymbol{\Theta}^{*} = \begin{bmatrix} -1.3 & -0.6 & 1 & 0\\ 0.4 & -1.7 & 0 & 1 \end{bmatrix}. \]

The full controller has eight scalar parameters. A diagonal-only feedback parameterization would force \( K_{x,12}=K_{x,21}=0 \) and could not reproduce the ideal matrix above. This simple check detects structural mismatch before any adaptive simulation.

10.1 Adaptation used in the software examples

Let \( \mathbf{e}=\mathbf{x}-\mathbf{x}_m \) and choose \( \mathbf{Q}=\mathbf{I}_2 \). For the diagonal reference matrix, the Lyapunov equation

\[ \mathbf{A}_m^{T}\mathbf{P} + \mathbf{P}\mathbf{A}_m = -\mathbf{I}_2 \]

has solution

\[ \mathbf{P} = \operatorname{diag} \left( \frac{1}{3},\frac{1}{4} \right). \]

The examples implement the normalized matrix updates

\[ \dot{\mathbf{K} }_x = -\frac{\gamma_x}{1+\mathbf{x}^{T}\mathbf{x} +\mathbf{r}^{T}\mathbf{r} } \mathbf{P}\mathbf{e}\mathbf{x}^{T}, \]

\[ \dot{\mathbf{K} }_r = -\frac{\gamma_r}{1+\mathbf{x}^{T}\mathbf{x} +\mathbf{r}^{T}\mathbf{r} } \mathbf{P}\mathbf{e}\mathbf{r}^{T}. \]

A Frobenius-norm projection bounds each gain matrix. The purpose is to demonstrate matrix parameterization, vectorization, logging, and implementation consistency. The exact assumptions required for a general MIMO MRAC stability theorem are addressed in Lesson 3.

Library coverage: NumPy provides dense matrix operations, Kronecker products, column-major vectorization, and CSV output. SciPy is not required because the selected diagonal Lyapunov equation has a closed-form solution.

Chapter23_Lesson2.py

"""
Chapter23_Lesson2.py
Parameterization of a 2x2 multivariable adaptive controller.

The plant is
    x_dot = A x + u
and the fixed reference model is
    x_m_dot = A_m x_m + r.

The controller is parameterized as
    u = K_x x + K_r r = Theta phi,
where Theta = [K_x K_r] and phi = [x^T r^T]^T.

This program:
1. verifies u = (phi^T kron I_m) vec(Theta);
2. simulates matrix-valued normalized gradient adaptation;
3. writes the trajectory to Chapter23_Lesson2_output.csv.

Dependencies:
    numpy
"""

from __future__ import annotations

import math
from pathlib import Path
import numpy as np


DT = 1.0e-3
FINAL_TIME = 20.0
GAMMA_X = 8.0
GAMMA_R = 8.0
GAIN_BOUND = 8.0


def project_frobenius(matrix: np.ndarray, bound: float) -> np.ndarray:
    """Project a matrix onto a Frobenius-norm ball."""
    norm = np.linalg.norm(matrix, ord="fro")
    if norm <= bound:
        return matrix
    return matrix * (bound / norm)


def reference(t: float) -> np.ndarray:
    """Bounded two-channel command with several frequencies."""
    return np.array(
        [
            math.sin(0.7 * t) + 0.4 * math.sin(1.3 * t),
            0.8 * math.cos(0.5 * t) + 0.25 * math.sin(1.1 * t),
        ],
        dtype=float,
    )


def verify_vectorization_identity() -> float:
    """Return the norm of the matrix-form/vector-form mismatch."""
    theta_matrix = np.array(
        [[1.2, -0.4, 0.7, 0.2], [-0.3, 0.9, -0.1, 1.1]],
        dtype=float,
    )
    phi = np.array([0.6, -0.8, 1.0, 0.25], dtype=float)

    u_matrix = theta_matrix @ phi
    # vec(.) is column-major in the standard Kronecker identity.
    theta_vector = theta_matrix.reshape(-1, order="F")
    psi = np.kron(phi.reshape(1, -1), np.eye(2))
    u_vector = psi @ theta_vector
    return float(np.linalg.norm(u_matrix - u_vector))


def simulate() -> tuple[np.ndarray, np.ndarray, np.ndarray, float]:
    """Simulate the plant, reference model, and adaptive matrix gains."""
    a = np.array([[-0.2, 0.6], [-0.4, -0.3]], dtype=float)
    a_m = np.array([[-1.5, 0.0], [0.0, -2.0]], dtype=float)

    # For Q = I, A_m^T P + P A_m = -Q has this diagonal solution.
    p = np.diag([1.0 / 3.0, 1.0 / 4.0])

    x = np.array([1.0, -0.5], dtype=float)
    x_m = np.zeros(2, dtype=float)
    k_x = np.zeros((2, 2), dtype=float)
    k_r = np.zeros((2, 2), dtype=float)

    steps = int(round(FINAL_TIME / DT))
    trajectory = np.zeros((steps + 1, 10), dtype=float)
    squared_error_integral = 0.0

    for step in range(steps + 1):
        t = step * DT
        r = reference(t)
        e = x - x_m
        u = k_x @ x + k_r @ r

        trajectory[step] = np.array(
            [t, x[0], x[1], x_m[0], x_m[1], e[0], e[1], u[0], u[1],
             np.linalg.norm(np.hstack((k_x, k_r)), ord="fro")]
        )

        if step == steps:
            break

        x_dot = a @ x + u
        x_m_dot = a_m @ x_m + r

        normalization = 1.0 + float(x @ x) + float(r @ r)
        direction = p @ e
        k_x_dot = -(GAMMA_X / normalization) * np.outer(direction, x)
        k_r_dot = -(GAMMA_R / normalization) * np.outer(direction, r)

        x = x + DT * x_dot
        x_m = x_m + DT * x_m_dot
        k_x = project_frobenius(k_x + DT * k_x_dot, GAIN_BOUND)
        k_r = project_frobenius(k_r + DT * k_r_dot, GAIN_BOUND)
        squared_error_integral += float(e @ e) * DT

    rms_error = math.sqrt(squared_error_integral / FINAL_TIME)
    return trajectory, k_x, k_r, rms_error


def main() -> None:
    mismatch = verify_vectorization_identity()
    trajectory, k_x, k_r, rms_error = simulate()

    output_path = Path(__file__).with_name("Chapter23_Lesson2_output.csv")
    header = "t,x1,x2,xm1,xm2,e1,e2,u1,u2,combined_gain_frobenius_norm"
    np.savetxt(output_path, trajectory, delimiter=",", header=header, comments="")

    np.set_printoptions(precision=6, suppress=True)
    print(f"Vectorization identity mismatch: {mismatch:.3e}")
    print(f"RMS tracking error: {rms_error:.6f}")
    print("Final K_x:")
    print(k_x)
    print("Final K_r:")
    print(k_r)
    print(f"Trajectory written to: {output_path}")


if __name__ == "__main__":
    main()

Library coverage: the implementation uses only the C++17 standard library and fixed-size arrays. For larger MIMO systems, Eigen is a common alternative for dynamic-size matrices, Kronecker products, decompositions, and linear solves.

Chapter23_Lesson2.cpp

/*
Chapter23_Lesson2.cpp
Fixed-size, dependency-free implementation of the same 2x2 adaptive
controller parameterization used in Chapter23_Lesson2.py.

Build:
    g++ -std=c++17 -O2 Chapter23_Lesson2.cpp -o Chapter23_Lesson2

Run:
    ./Chapter23_Lesson2
*/

#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>

using Vec2 = std::array<double, 2>;
using Mat2 = std::array<std::array<double, 2>, 2>;

constexpr double DT = 1.0e-3;
constexpr double FINAL_TIME = 20.0;
constexpr double GAMMA_X = 8.0;
constexpr double GAMMA_R = 8.0;
constexpr double GAIN_BOUND = 8.0;

Vec2 add(const Vec2& a, const Vec2& b) {
    return {a[0] + b[0], a[1] + b[1]};
}

Vec2 subtract(const Vec2& a, const Vec2& b) {
    return {a[0] - b[0], a[1] - b[1]};
}

Vec2 scale(const Vec2& a, double scalar) {
    return {scalar * a[0], scalar * a[1]};
}

Vec2 matVec(const Mat2& a, const Vec2& x) {
    return {
        a[0][0] * x[0] + a[0][1] * x[1],
        a[1][0] * x[0] + a[1][1] * x[1]
    };
}

double dot(const Vec2& a, const Vec2& b) {
    return a[0] * b[0] + a[1] * b[1];
}

Mat2 outer(const Vec2& a, const Vec2& b) {
    return { { {a[0] * b[0], a[0] * b[1]},
             {a[1] * b[0], a[1] * b[1]} } };
}

Mat2 addScaled(const Mat2& a, const Mat2& direction, double scalar) {
    Mat2 result = a;
    for (std::size_t i = 0; i < 2; ++i) {
        for (std::size_t j = 0; j < 2; ++j) {
            result[i][j] += scalar * direction[i][j];
        }
    }
    return result;
}

double frobeniusNorm(const Mat2& matrix) {
    double sum = 0.0;
    for (const auto& row : matrix) {
        for (double value : row) {
            sum += value * value;
        }
    }
    return std::sqrt(sum);
}

Mat2 projectFrobenius(const Mat2& matrix, double bound) {
    const double norm = frobeniusNorm(matrix);
    if (norm <= bound) {
        return matrix;
    }

    Mat2 projected = matrix;
    const double factor = bound / norm;
    for (auto& row : projected) {
        for (double& value : row) {
            value *= factor;
        }
    }
    return projected;
}

Vec2 reference(double t) {
    return {
        std::sin(0.7 * t) + 0.4 * std::sin(1.3 * t),
        0.8 * std::cos(0.5 * t) + 0.25 * std::sin(1.1 * t)
    };
}

double vectorizationMismatch() {
    // Theta is 2x4, stored by rows here.
    const std::array<std::array<double, 4>, 2> theta = { {
        { {1.2, -0.4, 0.7, 0.2} },
        { {-0.3, 0.9, -0.1, 1.1} }
    } };
    const std::array<double, 4> phi = { {0.6, -0.8, 1.0, 0.25} };

    Vec2 uMatrix = {0.0, 0.0};
    for (std::size_t i = 0; i < 2; ++i) {
        for (std::size_t j = 0; j < 4; ++j) {
            uMatrix[i] += theta[i][j] * phi[j];
        }
    }

    // Column-major vec(Theta) and Psi = phi^T kron I_2.
    std::array<double, 8> thetaVector{};
    for (std::size_t j = 0; j < 4; ++j) {
        for (std::size_t i = 0; i < 2; ++i) {
            thetaVector[2 * j + i] = theta[i][j];
        }
    }

    Vec2 uVector = {0.0, 0.0};
    for (std::size_t j = 0; j < 4; ++j) {
        uVector[0] += phi[j] * thetaVector[2 * j];
        uVector[1] += phi[j] * thetaVector[2 * j + 1];
    }

    const Vec2 difference = subtract(uMatrix, uVector);
    return std::sqrt(dot(difference, difference));
}

void printMatrix(const std::string& name, const Mat2& matrix) {
    std::cout << name << "\n";
    for (const auto& row : matrix) {
        std::cout << "  [" << std::setw(11) << row[0]
                  << ", " << std::setw(11) << row[1] << "]\n";
    }
}

int main() {
    try {
        const Mat2 a = { { {-0.2, 0.6}, {-0.4, -0.3} } };
        const Mat2 aM = { { {-1.5, 0.0}, {0.0, -2.0} } };
        const Mat2 p = { { {1.0 / 3.0, 0.0}, {0.0, 1.0 / 4.0} } };

        Vec2 x = {1.0, -0.5};
        Vec2 xM = {0.0, 0.0};
        Mat2 kX = { { {0.0, 0.0}, {0.0, 0.0} } };
        Mat2 kR = { { {0.0, 0.0}, {0.0, 0.0} } };

        std::ofstream csv("Chapter23_Lesson2_output.csv");
        if (!csv) {
            throw std::runtime_error("Cannot open Chapter23_Lesson2_output.csv");
        }
        csv << "t,x1,x2,xm1,xm2,e1,e2,u1,u2,combined_gain_frobenius_norm\n";
        csv << std::setprecision(12);

        const int steps = static_cast<int>(std::llround(FINAL_TIME / DT));
        double squaredErrorIntegral = 0.0;

        for (int step = 0; step <= steps; ++step) {
            const double t = step * DT;
            const Vec2 r = reference(t);
            const Vec2 e = subtract(x, xM);
            const Vec2 u = add(matVec(kX, x), matVec(kR, r));
            const double combinedNorm =
                std::sqrt(std::pow(frobeniusNorm(kX), 2.0) +
                          std::pow(frobeniusNorm(kR), 2.0));

            csv << t << ',' << x[0] << ',' << x[1] << ','
                << xM[0] << ',' << xM[1] << ','
                << e[0] << ',' << e[1] << ','
                << u[0] << ',' << u[1] << ','
                << combinedNorm << '\n';

            if (step == steps) {
                break;
            }

            const Vec2 xDot = add(matVec(a, x), u);
            const Vec2 xMDot = add(matVec(aM, xM), r);
            const double normalization = 1.0 + dot(x, x) + dot(r, r);
            const Vec2 direction = matVec(p, e);

            const Mat2 kXDot = outer(direction, x);
            const Mat2 kRDot = outer(direction, r);

            x = add(x, scale(xDot, DT));
            xM = add(xM, scale(xMDot, DT));
            kX = projectFrobenius(
                addScaled(kX, kXDot, -DT * GAMMA_X / normalization),
                GAIN_BOUND
            );
            kR = projectFrobenius(
                addScaled(kR, kRDot, -DT * GAMMA_R / normalization),
                GAIN_BOUND
            );
            squaredErrorIntegral += dot(e, e) * DT;
        }

        const double rmsError =
            std::sqrt(squaredErrorIntegral / FINAL_TIME);

        std::cout << std::fixed << std::setprecision(6);
        std::cout << "Vectorization identity mismatch: "
                  << vectorizationMismatch() << '\n';
        std::cout << "RMS tracking error: " << rmsError << '\n';
        printMatrix("Final K_x:", kX);
        printMatrix("Final K_r:", kR);
        std::cout << "Trajectory written to Chapter23_Lesson2_output.csv\n";
        return 0;
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
}

Library coverage: the example is dependency-free and uses primitive arrays. For general matrix dimensions, EJML or Apache Commons Math can replace the fixed-size helper methods.

Chapter23_Lesson2.java

/*
Chapter23_Lesson2.java
Dependency-free Java implementation of a 2x2 adaptive MIMO controller.

Compile:
    javac Chapter23_Lesson2.java

Run:
    java Chapter23_Lesson2
*/

import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.charset.StandardCharsets;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.Locale;

public final class Chapter23_Lesson2 {
    private static final double DT = 1.0e-3;
    private static final double FINAL_TIME = 20.0;
    private static final double GAMMA_X = 8.0;
    private static final double GAMMA_R = 8.0;
    private static final double GAIN_BOUND = 8.0;

    private Chapter23_Lesson2() {
    }

    private static double[] add(double[] a, double[] b) {
        return new double[] {a[0] + b[0], a[1] + b[1]};
    }

    private static double[] subtract(double[] a, double[] b) {
        return new double[] {a[0] - b[0], a[1] - b[1]};
    }

    private static double[] scale(double[] a, double scalar) {
        return new double[] {scalar * a[0], scalar * a[1]};
    }

    private static double[] matVec(double[][] a, double[] x) {
        return new double[] {
            a[0][0] * x[0] + a[0][1] * x[1],
            a[1][0] * x[0] + a[1][1] * x[1]
        };
    }

    private static double dot(double[] a, double[] b) {
        return a[0] * b[0] + a[1] * b[1];
    }

    private static double[][] outer(double[] a, double[] b) {
        return new double[][] {
            {a[0] * b[0], a[0] * b[1]},
            {a[1] * b[0], a[1] * b[1]}
        };
    }

    private static double[][] addScaled(
            double[][] matrix, double[][] direction, double scalar) {
        double[][] result = new double[2][2];
        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 2; ++j) {
                result[i][j] = matrix[i][j] + scalar * direction[i][j];
            }
        }
        return result;
    }

    private static double frobeniusNorm(double[][] matrix) {
        double sum = 0.0;
        for (double[] row : matrix) {
            for (double value : row) {
                sum += value * value;
            }
        }
        return Math.sqrt(sum);
    }

    private static double[][] projectFrobenius(
            double[][] matrix, double bound) {
        double norm = frobeniusNorm(matrix);
        if (norm <= bound) {
            return matrix;
        }

        double factor = bound / norm;
        double[][] result = new double[2][2];
        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 2; ++j) {
                result[i][j] = factor * matrix[i][j];
            }
        }
        return result;
    }

    private static double[] reference(double t) {
        return new double[] {
            Math.sin(0.7 * t) + 0.4 * Math.sin(1.3 * t),
            0.8 * Math.cos(0.5 * t) + 0.25 * Math.sin(1.1 * t)
        };
    }

    private static double vectorizationMismatch() {
        double[][] theta = {
            {1.2, -0.4, 0.7, 0.2},
            {-0.3, 0.9, -0.1, 1.1}
        };
        double[] phi = {0.6, -0.8, 1.0, 0.25};

        double[] uMatrix = {0.0, 0.0};
        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 4; ++j) {
                uMatrix[i] += theta[i][j] * phi[j];
            }
        }

        // Column-major vec(Theta).
        double[] thetaVector = new double[8];
        for (int j = 0; j < 4; ++j) {
            for (int i = 0; i < 2; ++i) {
                thetaVector[2 * j + i] = theta[i][j];
            }
        }

        double[] uVector = {0.0, 0.0};
        for (int j = 0; j < 4; ++j) {
            uVector[0] += phi[j] * thetaVector[2 * j];
            uVector[1] += phi[j] * thetaVector[2 * j + 1];
        }

        double[] difference = subtract(uMatrix, uVector);
        return Math.sqrt(dot(difference, difference));
    }

    private static void printMatrix(String name, double[][] matrix) {
        System.out.println(name);
        for (double[] row : matrix) {
            System.out.printf(Locale.US, "  [%11.6f, %11.6f]%n",
                    row[0], row[1]);
        }
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);

        double[][] a = { {-0.2, 0.6}, {-0.4, -0.3} };
        double[][] aM = { {-1.5, 0.0}, {0.0, -2.0} };
        double[][] p = { {1.0 / 3.0, 0.0}, {0.0, 1.0 / 4.0} };

        double[] x = {1.0, -0.5};
        double[] xM = {0.0, 0.0};
        double[][] kX = { {0.0, 0.0}, {0.0, 0.0} };
        double[][] kR = { {0.0, 0.0}, {0.0, 0.0} };

        int steps = (int) Math.round(FINAL_TIME / DT);
        double squaredErrorIntegral = 0.0;
        Path output = Path.of("Chapter23_Lesson2_output.csv");

        try (BufferedWriter writer = Files.newBufferedWriter(
                output, StandardCharsets.UTF_8)) {
            writer.write(
                "t,x1,x2,xm1,xm2,e1,e2,u1,u2,"
                + "combined_gain_frobenius_norm\n"
            );

            for (int step = 0; step <= steps; ++step) {
                double t = step * DT;
                double[] r = reference(t);
                double[] e = subtract(x, xM);
                double[] u = add(matVec(kX, x), matVec(kR, r));
                double combinedNorm = Math.sqrt(
                    Math.pow(frobeniusNorm(kX), 2.0)
                    + Math.pow(frobeniusNorm(kR), 2.0)
                );

                writer.write(String.format(
                    Locale.US,
                    "%.12g,%.12g,%.12g,%.12g,%.12g,"
                    + "%.12g,%.12g,%.12g,%.12g,%.12g%n",
                    t, x[0], x[1], xM[0], xM[1],
                    e[0], e[1], u[0], u[1], combinedNorm
                ));

                if (step == steps) {
                    break;
                }

                double[] xDot = add(matVec(a, x), u);
                double[] xMDot = add(matVec(aM, xM), r);
                double normalization = 1.0 + dot(x, x) + dot(r, r);
                double[] direction = matVec(p, e);

                double[][] kXDot = outer(direction, x);
                double[][] kRDot = outer(direction, r);

                x = add(x, scale(xDot, DT));
                xM = add(xM, scale(xMDot, DT));
                kX = projectFrobenius(
                    addScaled(
                        kX, kXDot, -DT * GAMMA_X / normalization
                    ),
                    GAIN_BOUND
                );
                kR = projectFrobenius(
                    addScaled(
                        kR, kRDot, -DT * GAMMA_R / normalization
                    ),
                    GAIN_BOUND
                );
                squaredErrorIntegral += dot(e, e) * DT;
            }
        } catch (IOException error) {
            System.err.println("I/O error: " + error.getMessage());
            System.exit(1);
        }

        double rmsError = Math.sqrt(
            squaredErrorIntegral / FINAL_TIME
        );
        System.out.printf(Locale.US,
                "Vectorization identity mismatch: %.3e%n",
                vectorizationMismatch());
        System.out.printf(Locale.US,
                "RMS tracking error: %.6f%n", rmsError);
        printMatrix("Final K_x:", kX);
        printMatrix("Final K_r:", kR);
        System.out.println("Trajectory written to: "
                + output.toAbsolutePath());
    }
}

MATLAB coverage: the script uses base MATLAB matrix syntax, kron, table output, and plotting. The Simulink script programmatically creates a compact fixed-step model containing a Clock source, a MATLAB Function block, and a workspace sink.

Chapter23_Lesson2.m

% Chapter23_Lesson2.m
% Parameterization and simulation of a 2x2 adaptive MIMO controller.
%
% The controller is u = Kx*x + Kr*r = Theta*phi.
% This script verifies the Kronecker-product vectorization identity,
% simulates normalized matrix-gain adaptation, and plots the result.
%
% Required software:
%   Base MATLAB
%
% Optional:
%   Simulink, using Chapter23_Lesson2_Simulink.m

clear; clc; close all;

dt = 1.0e-3;
finalTime = 20.0;
gammaX = 8.0;
gammaR = 8.0;
gainBound = 8.0;

A = [-0.2, 0.6; -0.4, -0.3];
Am = [-1.5, 0.0; 0.0, -2.0];
P = diag([1/3, 1/4]);

x = [1.0; -0.5];
xm = zeros(2, 1);
Kx = zeros(2, 2);
Kr = zeros(2, 2);

% Verify u = (phi' kron I_m) vec(Theta).
ThetaTest = [1.2, -0.4, 0.7, 0.2; ...
            -0.3, 0.9, -0.1, 1.1];
phiTest = [0.6; -0.8; 1.0; 0.25];
uMatrix = ThetaTest * phiTest;
uVector = kron(phiTest.', eye(2)) * ThetaTest(:);
vectorizationMismatch = norm(uMatrix - uVector);

steps = round(finalTime / dt);
logData = zeros(steps + 1, 10);
squaredErrorIntegral = 0.0;

for step = 0:steps
    t = step * dt;
    r = referenceSignal(t);
    e = x - xm;
    u = Kx * x + Kr * r;

    row = step + 1;
    logData(row, :) = [t, x.', xm.', e.', u.', ...
        norm([Kx, Kr], 'fro')];

    if step == steps
        break;
    end

    xDot = A * x + u;
    xmDot = Am * xm + r;

    normalization = 1.0 + x.' * x + r.' * r;
    direction = P * e;
    KxDot = -(gammaX / normalization) * direction * x.';
    KrDot = -(gammaR / normalization) * direction * r.';

    x = x + dt * xDot;
    xm = xm + dt * xmDot;
    Kx = projectFrobenius(Kx + dt * KxDot, gainBound);
    Kr = projectFrobenius(Kr + dt * KrDot, gainBound);
    squaredErrorIntegral = squaredErrorIntegral + e.' * e * dt;
end

rmsError = sqrt(squaredErrorIntegral / finalTime);

variableNames = {'t', 'x1', 'x2', 'xm1', 'xm2', ...
    'e1', 'e2', 'u1', 'u2', 'combined_gain_frobenius_norm'};
resultTable = array2table(logData, 'VariableNames', variableNames);
writetable(resultTable, 'Chapter23_Lesson2_output.csv');

fprintf('Vectorization identity mismatch: %.3e\n', ...
    vectorizationMismatch);
fprintf('RMS tracking error: %.6f\n', rmsError);
disp('Final Kx:');
disp(Kx);
disp('Final Kr:');
disp(Kr);

figure('Name', 'Chapter 23 Lesson 2: MIMO adaptive parameterization');
plot(logData(:, 1), logData(:, 2), 'LineWidth', 1.2);
hold on;
plot(logData(:, 1), logData(:, 4), '--', 'LineWidth', 1.2);
plot(logData(:, 1), logData(:, 3), 'LineWidth', 1.2);
plot(logData(:, 1), logData(:, 5), '--', 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('State / reference-model state');
legend('x_1', 'x_{m1}', 'x_2', 'x_{m2}', ...
    'Location', 'best');
title('Plant and reference-model trajectories');

figure('Name', 'Tracking error');
plot(logData(:, 1), logData(:, 6:7), 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('Tracking error');
legend('e_1', 'e_2', 'Location', 'best');
title('MIMO tracking error');

function r = referenceSignal(t)
    r = [sin(0.7 * t) + 0.4 * sin(1.3 * t); ...
         0.8 * cos(0.5 * t) + 0.25 * sin(1.1 * t)];
end

function matrix = projectFrobenius(matrix, bound)
    matrixNorm = norm(matrix, 'fro');
    if matrixNorm > bound
        matrix = matrix * (bound / matrixNorm);
    end
end

Chapter23_Lesson2_Simulink.m

% Chapter23_Lesson2_Simulink.m
% Programmatically build a compact Simulink realization of the
% Chapter 23, Lesson 2 adaptive MIMO example.
%
% Required software:
%   MATLAB and Simulink
%
% The model uses a Clock source, a MATLAB Function block containing the
% 2x2 plant/reference/adaptation recursion, and a To Workspace sink.

clear; clc;

if ~license('test', 'Simulink')
    error('Simulink is required to build this model.');
end

modelName = 'Chapter23_Lesson2_Model';

if bdIsLoaded(modelName)
    close_system(modelName, 0);
end

new_system(modelName);
open_system(modelName);

add_block('simulink/Sources/Clock', ...
    [modelName '/Clock'], ...
    'Position', [40, 90, 70, 120]);

add_block('simulink/User-Defined Functions/MATLAB Function', ...
    [modelName '/Adaptive MIMO Core'], ...
    'Position', [140, 55, 340, 155]);

add_block('simulink/Sinks/To Workspace', ...
    [modelName '/Trajectory'], ...
    'VariableName', 'Chapter23Lesson2SimulinkOutput', ...
    'SaveFormat', 'Array', ...
    'Position', [420, 85, 540, 125]);

add_line(modelName, 'Clock/1', 'Adaptive MIMO Core/1', ...
    'autorouting', 'on');
add_line(modelName, 'Adaptive MIMO Core/1', 'Trajectory/1', ...
    'autorouting', 'on');

root = sfroot;
chart = root.find('-isa', 'Stateflow.EMChart', ...
    'Path', [modelName '/Adaptive MIMO Core']);

chart.Script = sprintf([ ...
'function y = fcn(t)\n' ...
'%%#codegen\n' ...
'persistent x xm Kx Kr\n' ...
'if isempty(x)\n' ...
'    x = [1.0; -0.5];\n' ...
'    xm = zeros(2,1);\n' ...
'    Kx = zeros(2,2);\n' ...
'    Kr = zeros(2,2);\n' ...
'end\n' ...
'dt = 1.0e-3;\n' ...
'A = [-0.2, 0.6; -0.4, -0.3];\n' ...
'Am = [-1.5, 0.0; 0.0, -2.0];\n' ...
'P = diag([1/3, 1/4]);\n' ...
'gammaX = 8.0;\n' ...
'gammaR = 8.0;\n' ...
'gainBound = 8.0;\n' ...
'r = [sin(0.7*t) + 0.4*sin(1.3*t); ...\n' ...
'     0.8*cos(0.5*t) + 0.25*sin(1.1*t)];\n' ...
'e = x - xm;\n' ...
'u = Kx*x + Kr*r;\n' ...
'xDot = A*x + u;\n' ...
'xmDot = Am*xm + r;\n' ...
'normalization = 1.0 + x.''*x + r.''*r;\n' ...
'direction = P*e;\n' ...
'Kx = Kx - dt*(gammaX/normalization)*direction*x.'';\n' ...
'Kr = Kr - dt*(gammaR/normalization)*direction*r.'';\n' ...
'nx = norm(Kx, ''fro'');\n' ...
'if nx > gainBound\n' ...
'    Kx = Kx*(gainBound/nx);\n' ...
'end\n' ...
'nr = norm(Kr, ''fro'');\n' ...
'if nr > gainBound\n' ...
'    Kr = Kr*(gainBound/nr);\n' ...
'end\n' ...
'x = x + dt*xDot;\n' ...
'xm = xm + dt*xmDot;\n' ...
'y = [x; xm; e; u];\n' ...
'end\n' ...
]);

set_param(modelName, ...
    'StopTime', '20', ...
    'SolverType', 'Fixed-step', ...
    'Solver', 'FixedStepDiscrete', ...
    'FixedStep', '0.001');

save_system(modelName);
sim(modelName);

disp(['Model created and simulated: ', modelName, '.slx']);
disp(['Workspace variable: ', ...
    'Chapter23Lesson2SimulinkOutput']);

Library coverage: the notebook uses built-in KroneckerProduct, matrix operations, Reap and Sow for trajectory collection, CSV export, and ListLinePlot.

Chapter23_Lesson2.nb

Notebook[{
 Cell["Chapter 23, Lesson 2: Parameterization of Multivariable Controllers", "Title"],
 Cell["Matrix-gain and Kronecker-product implementation for a 2x2 adaptive MIMO example.", "Text"],
 Cell[BoxData[
"ClearAll[\"Global`*\"];
dt = 0.001; finalTime = 20.0; gammaX = 8.0; gammaR = 8.0; gainBound = 8.0;
a = { {-0.2, 0.6}, {-0.4, -0.3} }; am = { {-1.5, 0.0}, {0.0, -2.0} };
p = DiagonalMatrix[{1/3, 1/4}];
reference[t_] := {Sin[0.7 t] + 0.4 Sin[1.3 t],
  0.8 Cos[0.5 t] + 0.25 Sin[1.1 t]};
projectFrobenius[matrix_, bound_] := Module[{matrixNorm},
  matrixNorm = Norm[matrix, \"Frobenius\"];
  If[matrixNorm <= bound, matrix, matrix (bound/matrixNorm)]
];
thetaTest = { {1.2, -0.4, 0.7, 0.2}, {-0.3, 0.9, -0.1, 1.1} };
phiTest = {0.6, -0.8, 1.0, 0.25};
uMatrix = thetaTest . phiTest; thetaVector = Flatten[Transpose[thetaTest]];
psi = KroneckerProduct[{phiTest}, IdentityMatrix[2]];
uVector = Flatten[psi.thetaVector];
vectorizationMismatch = Norm[uMatrix - uVector];
x = {1.0, -0.5}; xm = {0.0, 0.0};
kx = ConstantArray[0.0, {2, 2}]; kr = ConstantArray[0.0, {2, 2}];
steps = Round[finalTime/dt]; squaredErrorIntegral = 0.0;
trajectory = Reap[
  Do[
    t = step dt; r = reference[t]; e = x - xm; u = kx.x + kr.r;
    Sow[Join[{t}, x, xm, e, u,
      {Norm[Join[kx, kr, 2], \"Frobenius\"]}]];
    If[step < steps,
      xDot = a.x + u; xmDot = am.xm + r;
      normalization = 1.0 + x.x + r.r; direction = p.e;
      kxDot = -(gammaX/normalization) Outer[Times, direction, x];
      krDot = -(gammaR/normalization) Outer[Times, direction, r];
      x = x + dt xDot; xm = xm + dt xmDot;
      kx = projectFrobenius[kx + dt kxDot, gainBound];
      kr = projectFrobenius[kr + dt krDot, gainBound];
      squaredErrorIntegral += e.e dt;
    ], {step, 0, steps}
  ]
][[2, 1]];
rmsError = Sqrt[squaredErrorIntegral/finalTime];
Export[\"Chapter23_Lesson2_output.csv\",
  Prepend[trajectory, {\"t\", \"x1\", \"x2\", \"xm1\", \"xm2\", \"e1\",
    \"e2\", \"u1\", \"u2\", \"combined_gain_frobenius_norm\"}]
];
Print[\"Vectorization identity mismatch: \", N[vectorizationMismatch]];
Print[\"RMS tracking error: \", N[rmsError]];
Print[\"Final Kx = \", MatrixForm[N[kx]]];
Print[\"Final Kr = \", MatrixForm[N[kr]]];
ListLinePlot[
  {trajectory[[All, {1, 2}]], trajectory[[All, {1, 4}]],
   trajectory[[All, {1, 3}]], trajectory[[All, {1, 5}]]},
  PlotLegends -> {\"x1\", \"xm1\", \"x2\", \"xm2\"},
  Frame -> True, FrameLabel -> {\"Time (s)\", \"State\"},
  PlotLabel -> \"Plant and Reference-Model Trajectories\"
]
"], "Input"]
}, WindowTitle -> "Chapter23_Lesson2", StyleDefinitions -> "Default.nb"
]        

Problem 1: Parameter Count and Kronecker Dimensions

A plant has \( n=5 \) states and \( m=3 \) inputs. The command vector has \( q=2 \) components. For \( \mathbf{u}=\mathbf{K}_x\mathbf{x} +\mathbf{K}_r\mathbf{r} \):

1. find the dimensions of both gain matrices;
2. find the number of scalar adaptive parameters;
3. find the dimensions of the Kronecker regressor.

Solution.

\[ \mathbf{K}_x\in\mathbb{R}^{3\times5}, \qquad \mathbf{K}_r\in\mathbb{R}^{3\times2}. \]

Therefore \( \boldsymbol{\Theta}\in\mathbb{R}^{3\times7} \) and contains

\[ 3(5+2)=21 \]

scalar parameters. Since

\[ \boldsymbol{\Psi} = \boldsymbol{\phi}^{T}\otimes\mathbf{I}_3, \]

its dimensions are

\[ \boldsymbol{\Psi}\in\mathbb{R}^{3\times21}. \]

Problem 2: Matching Feasibility for an Underactuated Plant

Let

\[ \mathbf{B} = \begin{bmatrix} 1\\ 0 \end{bmatrix}, \qquad \mathbf{A}_m-\mathbf{A} = \begin{bmatrix} -2 & 1\\ 0.5 & -1 \end{bmatrix}. \]

Determine whether a state-feedback gain \( \mathbf{K}_x\in\mathbb{R}^{1\times2} \) can satisfy exact matching.

Solution. The column space of \( \mathbf{B} \) is

\[ \operatorname{im}(\mathbf{B}) = \operatorname{span} \left\{ \begin{bmatrix} 1\\ 0 \end{bmatrix} \right\}. \]

Every column of \( \mathbf{B}\mathbf{K}_x \) must have zero second component. The desired correction has second row \( [0.5\ -1] \), so its columns do not lie in \( \operatorname{im}(\mathbf{B}) \). Exact matching is impossible with this input matrix and controller structure.

Problem 3: All Matching Gains for an Overactuated Plant

Let

\[ \mathbf{B} = \begin{bmatrix} 1 & 0 & 1\\ 0 & 1 & 1 \end{bmatrix}, \qquad \mathbf{D} = \begin{bmatrix} -1 & 0\\ 0 & -2 \end{bmatrix}. \]

Give the general form of all matrices \( \mathbf{K}\in\mathbb{R}^{3\times2} \) satisfying \( \mathbf{B}\mathbf{K}=\mathbf{D} \).

Solution. Since \( \mathbf{B} \) has full row rank,

\[ \mathbf{B}^{\dagger} = \mathbf{B}^{T} \left( \mathbf{B}\mathbf{B}^{T} \right)^{-1}. \]

All solutions are

\[ \mathbf{K} = \mathbf{B}^{\dagger}\mathbf{D} + \left( \mathbf{I}_3-\mathbf{B}^{\dagger}\mathbf{B} \right)\mathbf{Z}, \]

where \( \mathbf{Z}\in\mathbb{R}^{3\times2} \) is arbitrary. The second term allocates control effort in the null space of \( \mathbf{B} \) without changing the plant state derivative.

Problem 4: Reduced Symmetric Parameterization

For

\[ \mathbf{K}_x = \begin{bmatrix} \alpha_1 & \alpha_2\\ \alpha_2 & \alpha_1 \end{bmatrix}, \]

derive the regression \( \mathbf{u}=\boldsymbol{\Psi}_{\alpha} \boldsymbol{\alpha} \) for \( \mathbf{u}=\mathbf{K}_x\mathbf{x} \).

Solution.

\[ \mathbf{u} = \alpha_1 \begin{bmatrix} x_1\\ x_2 \end{bmatrix} + \alpha_2 \begin{bmatrix} x_2\\ x_1 \end{bmatrix}. \]

Therefore,

\[ \boldsymbol{\Psi}_{\alpha}(\mathbf{x}) = \begin{bmatrix} x_1 & x_2\\ x_2 & x_1 \end{bmatrix}, \qquad \boldsymbol{\alpha} = \begin{bmatrix} \alpha_1\\ \alpha_2 \end{bmatrix}. \]

The parameterization is identifiable at an instant only when \( x_1^2-x_2^2 \) is nonzero. Over time, the integrated information matrix may still become nonsingular even if the instantaneous matrix is singular at isolated times.

Problem 5: Coordinate Dependence of Diagonal Feedback

Let

\[ \mathbf{K}_x = \begin{bmatrix} k_1 & 0\\ 0 & k_2 \end{bmatrix}, \qquad \mathbf{x} = \begin{bmatrix} 1 & 1\\ 1 & -1 \end{bmatrix} \mathbf{z}. \]

Find \( \mathbf{K}_z \).

Solution.

\[ \mathbf{K}_z = \mathbf{K}_x\mathbf{T} = \begin{bmatrix} k_1 & k_1\\ k_2 & -k_2 \end{bmatrix}. \]

The transformed gain is generally dense. Hence, diagonal feedback is not an invariant concept under arbitrary state-coordinate changes.

Problem 6: Information-Matrix Rank

Suppose \( m=2 \) and the base regressor has four components. Over an experiment,

\[ \operatorname{rank} \left( \int_{0}^{T} \boldsymbol{\phi}\boldsymbol{\phi}^{T}\,dt \right) = 3. \]

What is the rank of the dense controller information matrix, and how many scalar parameter directions are unexcited?

Solution.

\[ \operatorname{rank}(\mathbf{G}_{\theta}) = m\times3 = 6. \]

The dense parameter vector has \( 2\times4=8 \) entries. Therefore, \( 8-6=2 \) scalar parameter directions are not distinguishable from the experiment.

A multivariable controller can be represented compactly as \( \mathbf{u}=\boldsymbol{\Theta}\boldsymbol{\phi} \) and equivalently as \( \mathbf{u} =(\boldsymbol{\phi}^{T}\otimes\mathbf{I}_m) \operatorname{vec}(\boldsymbol{\Theta}) \). Exact reference matching imposes algebraic range-space conditions that must be checked before adaptation. Dense matrices maximize flexibility, while masks, block structures, symmetry bases, and shared parameters reduce complexity and encode prior information. Fixed-filter dynamic regressors preserve linearity in adjustable parameters, whereas parameter-dependent controller states generally do not. Finally, parameter convergence depends on the information matrix of the chosen regressor and parameter basis; good tracking alone does not imply unique gain convergence.