In previous lessons, you learned how block-diagram algebra and signal flow graphs reduce interconnections into an equivalent transfer function. Real engineered systems, however, often contain:

• Multiple loops: nested feedback loops, interacting loops, and multiple summing junctions.
• Multiple signals: several inputs and outputs coupled through cross-channels.
• Directional dependencies: ordering matters (matrix multiplication is not commutative).

The correct language for this setting is the \( p \times m \) transfer-function matrix \( \mathbf{G}(s) \), mapping an input vector \( \mathbf{u}(s) \in \mathbb{R}^{m} \) to an output vector \( \mathbf{y}(s) \in \mathbb{R}^{p} \) by \( \mathbf{y}(s)=\mathbf{G}(s)\mathbf{u}(s) \).

flowchart LR
  R["r (m x 1)"] --> S1["Sum: e = r - H*y"]
  Y["y (p x 1)"] --> H["H(s) (m x p)"]
  H --> S1
  S1 --> K["K(s) (m x m)"]
  K --> G["G(s) (p x m)"]
  G --> Y
        

The diagram above depicts a standard MIMO negative-feedback interconnection. This lesson makes the diagram algebraic, proving when and how you may “collapse” the network into an equivalent transfer matrix.

Let \( \mathbf{u}(t)\in\mathbb{R}^{m} \) and \( \mathbf{y}(t)\in\mathbb{R}^{p} \) be input and output vectors of a linear time-invariant system (zero initial conditions). Their Laplace transforms satisfy

\[ \mathbf{y}(s)=\mathbf{G}(s)\mathbf{u}(s), \quad \mathbf{G}(s)\in\mathbb{R}(s)^{p\times m}. \]

The entry \( g_{ij}(s) \) is defined by the superposition experiment: set all input components to zero except the \( j \)th, and measure the \( i \)th output:

\[ g_{ij}(s) \;=\; \left.\frac{Y_i(s)}{U_j(s)}\right|_{U_k(s)=0 \text{ for } k\neq j}. \]

Dimension compatibility is enforced by matrix shapes: \( \mathbf{G}(s) \) must be \( p\times m \) so that the product \( \mathbf{G}(s)\mathbf{u}(s) \) is defined.

Remark (coupling): When \( \mathbf{G}(s) \) has nonzero off-diagonal terms, each input influences multiple outputs. This coupling is the central reason MIMO block diagram reduction must be done with matrix algebra rather than scalar “cancel-and-simplify” intuition.

This section proves the two simplest MIMO reduction identities. These proofs are important because they expose a key structural difference from SISO: matrix multiplication order matters.

3.1 Series (cascade) interconnection

Suppose \( \mathbf{G}_1(s) \) maps \( \mathbf{r}(s)\in\mathbb{R}^{k} \) to \( \mathbf{u}(s)\in\mathbb{R}^{m} \), and \( \mathbf{G}_2(s) \) maps \( \mathbf{u}(s)\in\mathbb{R}^{m} \) to \( \mathbf{y}(s)\in\mathbb{R}^{p} \):

\[ \mathbf{u}(s)=\mathbf{G}_1(s)\mathbf{r}(s),\qquad \mathbf{y}(s)=\mathbf{G}_2(s)\mathbf{u}(s). \]

Theorem (series rule): The equivalent transfer matrix is

\[ \mathbf{y}(s)=\mathbf{G}_{\mathrm{eq}}(s)\mathbf{r}(s),\qquad \mathbf{G}_{\mathrm{eq}}(s)=\mathbf{G}_2(s)\mathbf{G}_1(s). \]

Proof: Substitute \( \mathbf{u}(s)=\mathbf{G}_1(s)\mathbf{r}(s) \) into the second equation:

\[ \mathbf{y}(s)=\mathbf{G}_2(s)\mathbf{u}(s)=\mathbf{G}_2(s)\mathbf{G}_1(s)\mathbf{r}(s). \]

The product order is forced by signal flow. In general, \( \mathbf{G}_2\mathbf{G}_1\neq\mathbf{G}_1\mathbf{G}_2 \).

3.2 Parallel interconnection

Suppose two blocks share the same input \( \mathbf{r}(s)\in\mathbb{R}^{m} \) and their outputs add:

\[ \mathbf{y}(s)=\mathbf{G}_1(s)\mathbf{r}(s)+\mathbf{G}_2(s)\mathbf{r}(s). \]

Theorem (parallel rule):

\[ \mathbf{G}_{\mathrm{eq}}(s)=\mathbf{G}_1(s)+\mathbf{G}_2(s). \]

Proof: Factor \( \mathbf{r}(s) \) on the right:

\[ \mathbf{y}(s)=\big(\mathbf{G}_1(s)+\mathbf{G}_2(s)\big)\mathbf{r}(s). \]

The sum requires equal dimensions: both \( \mathbf{G}_1 \) and \( \mathbf{G}_2 \) must be the same size.

3.3 Non-commutativity as a design reality

Consider constant matrices (static gains) to isolate the linear-algebra point:

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

Then \( \mathbf{G}_2\mathbf{G}_1\neq\mathbf{G}_1\mathbf{G}_2 \), so swapping cascade order changes the system, even though both are “two blocks in series”.

We now derive the fundamental MIMO feedback identity. This identity is the matrix generalization of the familiar SISO formula \( \frac{G}{1+GH} \), but with a crucial structural detail: \( ( \mathbf{I}+\mathbf{G}\mathbf{H})^{-1} \) is a matrix inverse and depends on multiplication order.

4.1 Output feedback at the plant input

Let the plant be \( \mathbf{G}(s)\in\mathbb{R}(s)^{p\times m} \). Let the feedback block be \( \mathbf{H}(s)\in\mathbb{R}(s)^{m\times p} \). Consider the interconnection

\[ \mathbf{y}(s)=\mathbf{G}(s)\mathbf{u}(s),\qquad \mathbf{u}(s)=\mathbf{r}(s)-\mathbf{H}(s)\mathbf{y}(s). \]

Theorem (MIMO feedback reduction): If \( \mathbf{I}+\mathbf{G}(s)\mathbf{H}(s) \) is invertible, then

\[ \mathbf{y}(s)=\big(\mathbf{I}+\mathbf{G}(s)\mathbf{H}(s)\big)^{-1}\mathbf{G}(s)\,\mathbf{r}(s). \]

Proof: Substitute \( \mathbf{u}=\mathbf{r}-\mathbf{H}\mathbf{y} \) into \( \mathbf{y}=\mathbf{G}\mathbf{u} \):

\[ \mathbf{y} =\mathbf{G}\big(\mathbf{r}-\mathbf{H}\mathbf{y}\big) =\mathbf{G}\mathbf{r}-\mathbf{G}\mathbf{H}\mathbf{y}. \]

Bring terms in \( \mathbf{y} \) to the left:

\[ \big(\mathbf{I}+\mathbf{G}\mathbf{H}\big)\mathbf{y}=\mathbf{G}\mathbf{r}. \]

If \( \mathbf{I}+\mathbf{G}\mathbf{H} \) is invertible, multiply both sides by its inverse to obtain the result. Uniqueness follows from invertibility: if two solutions existed, their difference would lie in the nullspace of \( \mathbf{I}+\mathbf{G}\mathbf{H} \), which must be trivial.

4.2 Closed-loop, sensitivity, and complementary sensitivity matrices

Define the loop transfer matrix \( \mathbf{L}(s)=\mathbf{G}(s)\mathbf{H}(s) \). Then define:

\[ \mathbf{S}(s)=\big(\mathbf{I}+\mathbf{L}(s)\big)^{-1}, \qquad \mathbf{T}(s)=\mathbf{S}(s)\mathbf{G}(s). \]

Here \( \mathbf{S}(s) \) is the matrix that maps reference-equivalent signals to internal signals, and \( \mathbf{T}(s) \) maps \( \mathbf{r} \) to \( \mathbf{y} \) for the chosen interconnection. In MIMO, \( \mathbf{S} \) and \( \mathbf{T} \) are generally full matrices, even if \( \mathbf{G} \) or \( \mathbf{H} \) is sparse.

4.3 Controller-in-the-forward-path form

A very common form in control engineering is: controller \( \mathbf{K}(s)\in\mathbb{R}(s)^{m\times p} \) acts on an error \( \mathbf{e}(s)=\mathbf{r}(s)-\mathbf{y}(s) \) to produce \( \mathbf{u} \):

\[ \mathbf{u}(s)=\mathbf{K}(s)\big(\mathbf{r}(s)-\mathbf{y}(s)\big),\qquad \mathbf{y}(s)=\mathbf{G}(s)\mathbf{u}(s). \]

Substitution gives:

\[ \mathbf{y} =\mathbf{G}\mathbf{K}\big(\mathbf{r}-\mathbf{y}\big) =\mathbf{G}\mathbf{K}\mathbf{r}-\mathbf{G}\mathbf{K}\mathbf{y}, \quad \Rightarrow \quad \big(\mathbf{I}+\mathbf{G}\mathbf{K}\big)\mathbf{y}=\mathbf{G}\mathbf{K}\mathbf{r}. \]

Hence:

\[ \mathbf{y}=\big(\mathbf{I}+\mathbf{G}\mathbf{K}\big)^{-1}\mathbf{G}\mathbf{K}\,\mathbf{r}. \]

Note carefully the product ordering: \( \mathbf{G}\mathbf{K} \) arises because \( \mathbf{u}=\mathbf{K}(\cdot) \) is upstream of \( \mathbf{y}=\mathbf{G}\mathbf{u} \).

Multi-loop block diagrams can contain multiple internal signals. Instead of “guessing” reductions, we write a linear system in the Laplace domain and eliminate internal variables with matrix inversion. This method is conceptually aligned with signal-flow-graph reduction but uses linear algebra directly.

5.1 General internal-variable model

Collect all internal node signals into a vector \( \mathbf{z}(s)\in\mathbb{R}^{n} \). Suppose the block diagram implies:

\[ \mathbf{z}(s)=\mathbf{M}(s)\mathbf{z}(s)+\mathbf{N}(s)\mathbf{r}(s),\qquad \mathbf{y}(s)=\mathbf{P}(s)\mathbf{z}(s)+\mathbf{Q}(s)\mathbf{r}(s). \]

Theorem (elimination / well-posedness): If \( \mathbf{I}-\mathbf{M}(s) \) is invertible, then

\[ \mathbf{z}(s)=\big(\mathbf{I}-\mathbf{M}(s)\big)^{-1}\mathbf{N}(s)\mathbf{r}(s), \qquad \mathbf{y}(s)=\Big(\mathbf{P}(s)\big(\mathbf{I}-\mathbf{M}(s)\big)^{-1}\mathbf{N}(s)+\mathbf{Q}(s)\Big)\mathbf{r}(s). \]

Proof: Rearrange the first equation:

\[ \big(\mathbf{I}-\mathbf{M}\big)\mathbf{z}=\mathbf{N}\mathbf{r}. \]

If \( \mathbf{I}-\mathbf{M} \) is invertible, multiply by its inverse and substitute into \( \mathbf{y}=\mathbf{P}\mathbf{z}+\mathbf{Q}\mathbf{r} \).

Interpretation: the multi-loop network is “well-posed” (has a unique signal solution) exactly when \( \mathbf{I}-\mathbf{M}(s) \) is invertible for the relevant complex variable \( s \). In simulation tools, non-invertibility is closely related to algebraic-loop issues.

5.2 Example: two summing junctions (outer + inner loop)

Consider a nested-loop structure with two summing points: \( \mathbf{e}_1=\mathbf{r}-\mathbf{H}_1\mathbf{y} \), \( \mathbf{x}=\mathbf{K}_1\mathbf{e}_1 \), \( \mathbf{e}_2=\mathbf{x}-\mathbf{H}_2\mathbf{y} \), \( \mathbf{u}=\mathbf{K}_2\mathbf{e}_2 \), \( \mathbf{y}=\mathbf{G}\mathbf{u} \). The block diagram is:

flowchart LR
  R["r"] --> S1["Sum1: e1 = r - H1*y"]
  Y["y"] --> H1["H1(s)"] --> S1
  S1 --> K1["K1(s)"] --> S2["Sum2: e2 = x - H2*y"]
  Y --> H2["H2(s)"] --> S2
  S2 --> K2["K2(s)"] --> G["G(s)"] --> Y
        

Substitution gives one linear equation in \( \mathbf{y} \):

\[ \mathbf{y} =\mathbf{G}\mathbf{K}_2\big(\mathbf{K}_1(\mathbf{r}-\mathbf{H}_1\mathbf{y})-\mathbf{H}_2\mathbf{y}\big) =\mathbf{G}\mathbf{K}_2\mathbf{K}_1\mathbf{r} -\mathbf{G}\mathbf{K}_2\mathbf{K}_1\mathbf{H}_1\mathbf{y} -\mathbf{G}\mathbf{K}_2\mathbf{H}_2\mathbf{y}. \]

Move all \( \mathbf{y} \)-terms to the left:

\[ \Big(\mathbf{I}+\mathbf{G}\mathbf{K}_2\mathbf{K}_1\mathbf{H}_1+\mathbf{G}\mathbf{K}_2\mathbf{H}_2\Big)\mathbf{y} =\mathbf{G}\mathbf{K}_2\mathbf{K}_1\mathbf{r}. \]

If the bracketed matrix is invertible, the closed-loop transfer matrix from \( \mathbf{r} \) to \( \mathbf{y} \) follows immediately by inversion.

5.3 Practical note: frequency-by-frequency evaluation

Even if symbolic manipulation is hard, you can evaluate transfer matrices on the imaginary axis \( s=j\omega \) for \( \omega > 0 \). The closed-loop frequency response for a feedback loop is

\[ \mathbf{T}(j\omega)=\big(\mathbf{I}+\mathbf{G}(j\omega)\mathbf{H}(j\omega)\big)^{-1}\mathbf{G}(j\omega), \qquad \omega > 0. \]

This is exactly what the language labs implement: evaluate each transfer-function entry at \( j\omega \), assemble matrices, invert numerically, and compute \( \mathbf{T}(j\omega) \).

A signal flow graph (SFG) is built from scalar node variables. For MIMO, the “inputs” are the components \( R_1,\dots,R_m \) and the “outputs” are \( Y_1,\dots,Y_p \). The transfer matrix entry \( g_{ij}(s) \) is computed by exciting only \( R_j \) and setting the other inputs to zero, then applying Mason’s Gain Formula to find \( Y_i/R_j \).

Proposition (MIMO via superposition): For an LTI network represented by an SFG, the transfer matrix satisfies

\[ \mathbf{G}(s)=\begin{bmatrix} \dfrac{Y_1}{R_1} & \cdots & \dfrac{Y_1}{R_m}\\ \vdots & \ddots & \vdots\\ \dfrac{Y_p}{R_1} & \cdots & \dfrac{Y_p}{R_m} \end{bmatrix}, \]

where each ratio is computed on the same SFG but with all inputs except the selected one set to zero.

Proof sketch (sufficient for this course stage):

1. Linearity implies superposition: the response to \( \mathbf{r}=\sum_{j=1}^m \mathbf{e}_j R_j \) is the sum of responses to each single-input excitation.
2. By definition, when only \( R_j \) is applied, the output component satisfies \( Y_i=g_{ij}R_j \).
3. Therefore, for a general input vector, \( Y_i=\sum_{j=1}^m g_{ij}R_j \), which is exactly \( \mathbf{y}=\mathbf{G}\mathbf{r} \).

Computationally, Mason’s formula is performed per input-output pair; the matrix assembly is bookkeeping.

This lab computes the closed-loop frequency response \( \mathbf{T}(j\omega)=\big(\mathbf{I}+\mathbf{G}(j\omega)\mathbf{H}(j\omega)\big)^{-1}\mathbf{G}(j\omega) \) using only complex arithmetic and matrix inversion. This avoids needing symbolic rational matrix algebra.

Recommended libraries: numpy for matrices and scipy for polynomial utilities. Many courses also use python-control (installable via pip install control) for direct MIMO feedback, but the implementation below is “from scratch” at the frequency-response level.

import numpy as np

def eval_poly(coeffs, s):
    """
    Evaluate polynomial c[0]*s^n + ... + c[n] at complex s.
    coeffs: list or 1D array, highest power first.
    """
    y = 0j
    for c in coeffs:
        y = y * s + c
    return y

def eval_tf(num, den, s):
    """Evaluate transfer function num(s)/den(s) at complex s."""
    return eval_poly(num, s) / eval_poly(den, s)

# Example 2x2 plant G(s): each entry is a proper transfer function
# G11 = 1/(s+1), G12 = 2/(s+2), G21 = 1/(s+3), G22 = 1/(s+1)
G_num = [
    [[1.0], [2.0]],   # row 1 numerators
    [[1.0], [1.0]]    # row 2 numerators
]
G_den = [
    [[1.0, 1.0], [1.0, 2.0]],  # row 1 denominators
    [[1.0, 3.0], [1.0, 1.0]]   # row 2 denominators
]

# Unity output feedback: H(s) = I (dimensions m x p = 2 x 2)
# So L(jw) = G(jw)*I = G(jw), and T(jw) = (I + G(jw))^{-1} * G(jw)
def G_eval(s):
    G = np.zeros((2, 2), dtype=complex)
    for i in range(2):
        for j in range(2):
            G[i, j] = eval_tf(G_num[i][j], G_den[i][j], s)
    return G

def T_eval(jw):
    s = 1j * jw
    G = G_eval(s)
    I = np.eye(2, dtype=complex)
    T = np.linalg.inv(I + G) @ G
    return T

# Evaluate on a frequency grid
w_grid = np.logspace(-2, 2, 9)  # rad/s
for w in w_grid:
    T = T_eval(w)
    # Print magnitudes of the 4 closed-loop channels
    mags = np.abs(T)
    print(f"w={w:8.4f}  |T|=\n{mags}\n")
      

What to observe: even with unity feedback, the off-diagonal elements of \( \mathbf{T}(j\omega) \) typically remain nonzero, reflecting closed-loop coupling.

In C++, we implement the same frequency-response computation using:

• std::complex<double> for complex arithmetic,
• Eigen for complex matrix inversion and multiplication.

#include <complex>
#include <iostream>
#include <vector>
#include <Eigen/Dense>

using cd = std::complex<double>;

cd eval_poly(const std::vector<double>& c, const cd& s) {
    // Horner: c[0]*s^n + ... + c[n]
    cd y = cd(0.0, 0.0);
    for (double a : c) {
        y = y * s + a;
    }
    return y;
}

cd eval_tf(const std::vector<double>& num,
           const std::vector<double>& den,
           const cd& s) {
    return eval_poly(num, s) / eval_poly(den, s);
}

int main() {
    // Example 2x2 plant:
    // G11 = 1/(s+1), G12 = 2/(s+2), G21 = 1/(s+3), G22 = 1/(s+1)
    std::vector<double> n11{1.0}, d11{1.0, 1.0};
    std::vector<double> n12{2.0}, d12{1.0, 2.0};
    std::vector<double> n21{1.0}, d21{1.0, 3.0};
    std::vector<double> n22{1.0}, d22{1.0, 1.0};

    std::vector<double> w_grid{0.01, 0.1, 1.0, 10.0, 100.0};

    for (double w : w_grid) {
        cd s(0.0, w); // s = j*w

        Eigen::Matrix<cd, 2, 2> G;
        G(0,0) = eval_tf(n11, d11, s);
        G(0,1) = eval_tf(n12, d12, s);
        G(1,0) = eval_tf(n21, d21, s);
        G(1,1) = eval_tf(n22, d22, s);

        Eigen::Matrix<cd, 2, 2> I = Eigen::Matrix<cd, 2, 2>::Identity();

        // Unity feedback: T = (I + G)^(-1) * G
        Eigen::Matrix<cd, 2, 2> T = (I + G).inverse() * G;

        std::cout << "w=" << w << "\n";
        std::cout << "|T|=\n";
        std::cout << std::abs(T(0,0)) << " " << std::abs(T(0,1)) << "\n";
        std::cout << std::abs(T(1,0)) << " " << std::abs(T(1,1)) << "\n\n";
    }
    return 0;
}
      

This computation mirrors the mathematical identity and is a robust way to validate multi-loop reductions: compute both the “reduced” and “unreduced” forms numerically and compare \( \mathbf{T}(j\omega) \).

In Java, you can implement the same approach using:

• org.apache.commons.math3.complex.Complex for complex arithmetic,
• a matrix library (e.g., EJML complex matrices) for inversion.

Below is an educational implementation using plain 2x2 complex matrix inversion (from scratch) to avoid relying on a specific complex-matrix package. This keeps the core idea transparent.

import org.apache.commons.math3.complex.Complex;

public class MimoFeedbackFRF {

    static Complex evalPoly(double[] c, Complex s) {
        Complex y = Complex.ZERO;
        for (double a : c) {
            y = y.multiply(s).add(a);
        }
        return y;
    }

    static Complex evalTF(double[] num, double[] den, Complex s) {
        return evalPoly(num, s).divide(evalPoly(den, s));
    }

    // 2x2 matrix utilities (Complex)
    static Complex[][] add(Complex[][] A, Complex[][] B) {
        Complex[][] C = new Complex[2][2];
        for (int i=0;i<2;i++) for (int j=0;j<2;j++) C[i][j] = A[i][j].add(B[i][j]);
        return C;
    }

    static Complex[][] mul(Complex[][] A, Complex[][] B) {
        Complex[][] C = new Complex[2][2];
        for (int i=0;i<2;i++) {
            for (int j=0;j<2;j++) {
                C[i][j] = Complex.ZERO;
                for (int k=0;k<2;k++) C[i][j] = C[i][j].add(A[i][k].multiply(B[k][j]));
            }
        }
        return C;
    }

    static Complex[][] inv2x2(Complex[][] A) {
        // A^{-1} = (1/det) * [[d, -b],[-c, a]]
        Complex a = A[0][0], b = A[0][1], c = A[1][0], d = A[1][1];
        Complex det = a.multiply(d).subtract(b.multiply(c));
        Complex[][] B = new Complex[2][2];
        B[0][0] = d.divide(det);
        B[0][1] = b.negate().divide(det);
        B[1][0] = c.negate().divide(det);
        B[1][1] = a.divide(det);
        return B;
    }

    static double abs(Complex z) { return z.abs(); }

    public static void main(String[] args) {
        // G11 = 1/(s+1), G12 = 2/(s+2), G21 = 1/(s+3), G22 = 1/(s+1)
        double[] n11 = {1.0}, d11 = {1.0, 1.0};
        double[] n12 = {2.0}, d12 = {1.0, 2.0};
        double[] n21 = {1.0}, d21 = {1.0, 3.0};
        double[] n22 = {1.0}, d22 = {1.0, 1.0};

        double[] wGrid = {0.01, 0.1, 1.0, 10.0, 100.0};

        for (double w : wGrid) {
            Complex s = new Complex(0.0, w); // j*w

            Complex[][] G = new Complex[2][2];
            G[0][0] = evalTF(n11, d11, s);
            G[0][1] = evalTF(n12, d12, s);
            G[1][0] = evalTF(n21, d21, s);
            G[1][1] = evalTF(n22, d22, s);

            Complex[][] I = {
                {Complex.ONE, Complex.ZERO},
                {Complex.ZERO, Complex.ONE}
            };

            // Unity feedback: T = (I + G)^{-1} * G
            Complex[][] IpG = add(I, G);
            Complex[][] T = mul(inv2x2(IpG), G);

            System.out.println("w=" + w);
            System.out.println("|T|=");
            System.out.println(abs(T[0][0]) + " " + abs(T[0][1]));
            System.out.println(abs(T[1][0]) + " " + abs(T[1][1]));
            System.out.println();
        }
    }
}
      

MATLAB’s Control System Toolbox supports transfer-function matrices directly. You can compute MIMO feedback using feedback and multi-loop connections using connect, sumblk, and named I/O.

% MIMO transfer-function matrix example
s = tf('s');

G11 = 1/(s+1);
G12 = 2/(s+2);
G21 = 1/(s+3);
G22 = 1/(s+1);
G = [G11 G12; G21 G22];

% Unity feedback (H = I): T = (I + G)^(-1) * G  via feedback(G, I)
I2 = eye(2);
T = feedback(G, I2);   % negative feedback by default

% Frequency response check
w = logspace(-2, 2, 9);
Tjw = freqresp(T, w);  % size: 2x2xN
mag11 = squeeze(abs(Tjw(1,1,:)));

disp('Magnitudes of T11(jw):');
disp([w(:) mag11]);

% Nested-loop style example from Section 5:
% e1 = r - H1*y; x = K1*e1; e2 = x - H2*y; u = K2*e2; y = G*u
K1 = diag([5, 3]);            % example static controller gains
K2 = diag([2, 1]);            % inner controller gains
H1 = [1 0; 0 1];              % outer measurement matrix
H2 = [0.5 0; 0 0.2];          % inner measurement matrix

% Use algebraic reduction by building the combined equation:
% (I + G*K2*K1*H1 + G*K2*H2) y = G*K2*K1 r
I2 = eye(2);
T_nested = inv(I2 + G*K2*K1*H1 + G*K2*H2) * (G*K2*K1);

% Compare frequency response numerically
Tn_jw = freqresp(T_nested, w);
disp('Magnitudes of T_nested(1,2)(jw):');
disp([w(:) squeeze(abs(Tn_jw(1,2,:)))]);
      

Simulink note: For MIMO transfer functions, a convenient approach is to define G in the workspace and use an LTI System block (which accepts MIMO LTI objects), together with Sum and Gain blocks for matrix gains K, H. Keep signal dimensions explicit (vector widths) to avoid unintended scalar expansion.

Mathematica supports MIMO transfer-function models and feedback composition using system-model interconnection functions.

(* Define Laplace variable s *)
s = Global`s;

(* 2x2 transfer-function matrix G(s) *)
G11 = 1/(s + 1);
G12 = 2/(s + 2);
G21 = 1/(s + 3);
G22 = 1/(s + 1);

G = TransferFunctionModel[{ {G11, G12}, {G21, G22} }, s];

(* Unity feedback: H = I *)
H = TransferFunctionModel[IdentityMatrix[2], s];

(* Closed-loop: T = feedback(G, H) *)
T = SystemsModelFeedbackConnect[G, H];

(* Frequency response at selected frequencies *)
wGrid = {0.01, 0.1, 1, 10, 100};
Tjw = Table[
  {w, FrequencyResponse[T, I*w]},
  {w, wGrid}
];

Tjw
      

The following problems reinforce MIMO interconnection algebra, invertibility (well-posedness), and multi-loop reduction.

Problem 1 (MIMO feedback reduction): Given \( \mathbf{y}=\mathbf{G}\mathbf{u} \) and \( \mathbf{u}=\mathbf{r}-\mathbf{H}\mathbf{y} \), derive \( \mathbf{y} \) in terms of \( \mathbf{r} \), and state the condition for uniqueness.

Solution: Substitute and rearrange:

\[ \mathbf{y}=\mathbf{G}(\mathbf{r}-\mathbf{H}\mathbf{y}) =\mathbf{G}\mathbf{r}-\mathbf{G}\mathbf{H}\mathbf{y} \;\Rightarrow\; (\mathbf{I}+\mathbf{G}\mathbf{H})\mathbf{y}=\mathbf{G}\mathbf{r}. \]

If \( \mathbf{I}+\mathbf{G}\mathbf{H} \) is invertible, then

\[ \mathbf{y}=(\mathbf{I}+\mathbf{G}\mathbf{H})^{-1}\mathbf{G}\mathbf{r}. \]

Uniqueness holds exactly when \( \mathbf{I}+\mathbf{G}\mathbf{H} \) is invertible.

Problem 2 (static MIMO closed-loop gain): Let \( \mathbf{G}_0=\begin{bmatrix}2&1\\0&1\end{bmatrix} \) and \( \mathbf{H}_0=\begin{bmatrix}1&0\\1&1\end{bmatrix} \). Compute \( \mathbf{T}_0=(\mathbf{I}+\mathbf{G}_0\mathbf{H}_0)^{-1}\mathbf{G}_0 \).

Solution: First compute

\[ \mathbf{I}+\mathbf{G}_0\mathbf{H}_0 = \begin{bmatrix}1&0\\0&1\end{bmatrix} + \begin{bmatrix}2&1\\0&1\end{bmatrix} \begin{bmatrix}1&0\\1&1\end{bmatrix} = \begin{bmatrix}4&1\\1&2\end{bmatrix}. \]

Invert the 2x2 matrix using \( \begin{bmatrix}a&b\\c&d\end{bmatrix}^{-1}=\frac{1}{ad-bc}\begin{bmatrix}d&-b\\-c&a\end{bmatrix} \). Here \( ad-bc=4\cdot 2-1\cdot 1=7 \), so

\[ (\mathbf{I}+\mathbf{G}_0\mathbf{H}_0)^{-1} = \frac{1}{7}\begin{bmatrix}2&-1\\-1&4\end{bmatrix}. \]

Multiply by \( \mathbf{G}_0 \):

\[ \mathbf{T}_0 = \frac{1}{7}\begin{bmatrix}2&-1\\-1&4\end{bmatrix} \begin{bmatrix}2&1\\0&1\end{bmatrix} = \begin{bmatrix}\frac{4}{7}&\frac{1}{7}\\-\frac{2}{7}&\frac{3}{7}\end{bmatrix}. \]

Problem 3 (non-commutativity in cascade): Using \( \mathbf{G}_1=\begin{bmatrix}1&1\\0&1\end{bmatrix} \) and \( \mathbf{G}_2=\begin{bmatrix}1&0\\1&1\end{bmatrix} \), compute \( \mathbf{G}_2\mathbf{G}_1 \) and \( \mathbf{G}_1\mathbf{G}_2 \) and compare.

Solution:

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

Since the products differ, cascade order changes the system even for static gains.

Problem 4 (2x2 symbolic closed-loop for unity feedback): Let \( \mathbf{G}=\begin{bmatrix}g_{11}&g_{12}\\g_{21}&g_{22}\end{bmatrix} \) and unity feedback \( \mathbf{H}=\mathbf{I} \). Derive an explicit formula for \( \mathbf{T}=(\mathbf{I}+\mathbf{G})^{-1}\mathbf{G} \).

Solution: Compute

\[ \mathbf{I}+\mathbf{G} = \begin{bmatrix}1+g_{11}&g_{12}\\g_{21}&1+g_{22}\end{bmatrix}, \quad \Delta=(1+g_{11})(1+g_{22})-g_{12}g_{21}. \]

Assuming \( \Delta\neq 0 \), the inverse is

\[ (\mathbf{I}+\mathbf{G})^{-1} =\frac{1}{\Delta}\begin{bmatrix}1+g_{22}&-g_{12}\\-g_{21}&1+g_{11}\end{bmatrix}. \]

Thus

\[ \mathbf{T} =\frac{1}{\Delta} \begin{bmatrix}1+g_{22}&-g_{12}\\-g_{21}&1+g_{11}\end{bmatrix} \begin{bmatrix}g_{11}&g_{12}\\g_{21}&g_{22}\end{bmatrix}. \]

Multiplying gives each closed-loop channel explicitly; for example,

\[ T_{11}=\frac{(1+g_{22})g_{11}-g_{12}g_{21}}{\Delta}, \qquad T_{12}=\frac{(1+g_{22})g_{12}-g_{12}g_{22}}{\Delta}=\frac{g_{12}}{\Delta}. \]

This shows how coupling terms \( g_{12}g_{21} \) enter both numerator and denominator.

Problem 5 (well-posedness of a multi-loop internal system): Consider the internal-variable form \( \mathbf{z}=\mathbf{M}\mathbf{z}+\mathbf{N}\mathbf{r} \), \( \mathbf{y}=\mathbf{P}\mathbf{z} \). Prove that the network has a unique internal solution iff \( \mathbf{I}-\mathbf{M} \) is invertible.

Solution:

If \( \mathbf{I}-\mathbf{M} \) is invertible, then \( (\mathbf{I}-\mathbf{M})\mathbf{z}=\mathbf{N}\mathbf{r} \) implies the unique solution \( \mathbf{z}=(\mathbf{I}-\mathbf{M})^{-1}\mathbf{N}\mathbf{r} \).

Conversely, if \( \mathbf{I}-\mathbf{M} \) is not invertible, then there exists a nonzero \( \mathbf{v}\neq \mathbf{0} \) such that \( (\mathbf{I}-\mathbf{M})\mathbf{v}=\mathbf{0} \). For any particular solution \( \mathbf{z}_0 \), \( \mathbf{z}_0+\alpha \mathbf{v} \) is also a solution (for all scalars \( \alpha \)), so uniqueness fails.

We introduced transfer-function matrices as the correct representation for MIMO systems in block diagrams and signal flow graphs. Using rigorous linear-algebra arguments, we proved MIMO series and parallel rules, derived the single-loop feedback reduction \( (\mathbf{I}+\mathbf{G}\mathbf{H})^{-1}\mathbf{G} \), and developed a general internal-variable elimination method for multi-loop networks. The programming labs demonstrated a practical workflow: evaluate each transfer entry at \( s=j\omega \), assemble matrices, and compute closed-loop frequency responses by inversion.