The preceding chapters developed a complete finite-dimensional state-space vocabulary: states, inputs, outputs, modes, controllability, observability, realization, state feedback, integral action, and quadratic performance. Advanced control courses do not replace this vocabulary; they add optimization, uncertainty, constraints, multivariable interaction measures, or estimation theory on top of it.

We use the continuous-time LTI model \( \dot{\mathbf{x} }=A\mathbf{x}+B\mathbf{u} \), \( \mathbf{y}=C\mathbf{x}+D\mathbf{u} \) as the reference form throughout the lesson. A state-feedback controller has \( \mathbf{u}=-K\mathbf{x} \), giving the closed-loop state matrix \( A_{\mathrm{cl} }=A-BK \).

\[ \boxed{ \begin{aligned} &\text{Modern Control supplies the state-space language; advanced control courses}\\ &\text{add objectives, uncertainty, constraints, and estimation architecture.} \end{aligned} } \]

flowchart LR
  A["Modern Control core"] --> B["Optimal control: \nchoose K by minimizing cost"]
  A --> C["Robust control: \nguarantee stability \nunder uncertainty"]
  A --> D["MPC: optimize over a future \nhorizon with constraints"]
  A --> E["Multivariable control: \nmanage input-output coupling"]
  A --> F["Observer design: \nestimate missing states \nfrom measurements"]
  B --> G["Advanced state-space design"]
  C --> G
  D --> G
  E --> G
  F --> G
        

Pole placement asks for a feedback matrix \(K\) such that the eigenvalues of \(A-BK\) equal prescribed locations. Optimal control instead specifies a performance index and derives \(K\) from that index. The infinite-horizon LQR problem minimizes

\[ J(\mathbf{x}_0,\mathbf{u})=\int_0^\infty \left(\mathbf{x}(t)^\top Q\mathbf{x}(t)+ \mathbf{u}(t)^\top R\mathbf{u}(t)\right)\,dt, \quad Q\succeq 0,\quad R\succ 0. \]

If \((A,B)\) is stabilizable and \((Q^{1/2},A)\) is detectable, the stabilizing solution \(P=P^\top\succeq 0\) of the algebraic Riccati equation satisfies

\[ A^\top P+PA-PBR^{-1}B^\top P+Q=0,\qquad K=R^{-1}B^\top P. \]

The closed-loop poles are not chosen directly. They emerge from \(Q\) and \(R\). Large state weights increase pressure to reduce corresponding states; large input weights penalize actuator effort.

Stability proof by a Lyapunov argument. Define

\[ V(\mathbf{x})=\mathbf{x}^\top P\mathbf{x}. \]

Under \(\mathbf{u}=-K\mathbf{x}\), \(\dot{\mathbf{x} }=(A-BK)\mathbf{x}\). Therefore,

\[ \begin{aligned} \dot{V} &=\mathbf{x}^\top\left((A-BK)^\top P+P(A-BK)\right)\mathbf{x} \\ &= -\mathbf{x}^\top(Q+K^\top RK)\mathbf{x}. \end{aligned} \]

If the unpenalized modes are detectable, then \(\dot{V}<0\) along all nonzero closed-loop trajectories that matter for stability, so the LQR law is stabilizing. This proof shows why LQR is a systematic continuation of Lyapunov methods rather than a separate design philosophy.

Earlier lessons treated the nominal model \((A,B,C,D)\) as exact. Robust control begins when the model is only approximate. A simple additive state-matrix uncertainty is

\[ \dot{\mathbf{x} }=(A+\Delta_A)\mathbf{x}+B\mathbf{u},\qquad \mathbf{u}=-K\mathbf{x}. \]

The uncertain closed-loop matrix is \(A_{\Delta}=A-BK+\Delta_A\). A Lyapunov sufficient condition for robust stability over an uncertainty set \(\mathcal{D}\) is the existence of one \(P=P^\top\succ 0\) such that

\[ A_{\Delta}^{\top}P+PA_{\Delta}\prec 0,\qquad \forall\,\Delta_A\in\mathcal{D}. \]

This is the common-quadratic-Lyapunov viewpoint. In frequency-domain robust control, the same idea becomes a gain condition. For a stable nominal closed loop perturbed by a stable uncertainty block \(\Delta\), a small-gain sufficient condition has the form

\[ \left\|G_{\mathrm{cl} }(s)\Delta(s)\right\|_{\infty}<1. \]

Thus robust control extends Modern Control by asking not only whether \(A-BK\) is stable, but whether stability survives admissible modeling errors, parameter variation, neglected dynamics, and actuator/sensor perturbations.

MPC uses the state-space model as a prediction model. For the discrete system \(\mathbf{x}_{k+1}=A_d\mathbf{x}_k+B_d\mathbf{u}_k\), the finite-horizon quadratic problem is

\[ \begin{aligned} \min_{\mathbf{u}_0,\ldots,\mathbf{u}_{N-1} }\quad &\sum_{k=0}^{N-1} \left(\mathbf{x}_k^\top Q\mathbf{x}_k+ \mathbf{u}_k^\top R\mathbf{u}_k\right) +\mathbf{x}_N^\top Q_f\mathbf{x}_N \\ \text{subject to}\quad &\mathbf{x}_{k+1}=A_d\mathbf{x}_k+B_d\mathbf{u}_k, \\ &\mathbf{x}_k\in\mathcal{X},\qquad \mathbf{u}_k\in\mathcal{U}. \end{aligned} \]

MPC applies only the first computed input and then resolves the problem at the next sample. This is why MPC is also called receding-horizon control. Without constraints, the finite horizon problem reduces to a Riccati recursion:

\[ \begin{aligned} P_N&=Q_f,\\ K_k&=(R+B_d^\top P_{k+1}B_d)^{-1}B_d^\top P_{k+1}A_d,\\ P_k&=Q+A_d^\top P_{k+1}(A_d-B_dK_k). \end{aligned} \]

The constrained case is the distinguishing feature of MPC. State and input limits turn the problem into a quadratic program. Therefore, modern MPC is a direct extension of state-space simulation, controllability, quadratic performance, and numerical optimization.

In a MIMO system, inputs and outputs interact through matrices rather than scalar transfer functions. The transfer matrix is

\[ G(s)=C(sI-A)^{-1}B+D. \]

Pole placement and LQR already operate naturally on MIMO systems: \(B\) may have several columns and \(K\) may contain several rows. However, multivariable control adds questions about input-output directions, singular values, cross-coupling, decentralized structure, and robustness of directions.

\[ \bar{\sigma}(G(j\omega))= \sqrt{\lambda_{\max}(G(j\omega)^*G(j\omega))},\qquad \underline{\sigma}(G(j\omega))= \sqrt{\lambda_{\min}(G(j\omega)^*G(j\omega))}. \]

The largest singular value describes the greatest amplification from an input direction to an output direction at frequency \(\omega\). The smallest singular value describes the weakest direction. Modern Control prepares students for this by treating state, input, and output spaces as vector spaces rather than scalar signal paths.

\[ \boxed{ \begin{aligned} &\text{SISO design moves poles; MIMO design also shapes directions, }\\ &\text{gains, interactions, and allocation of control authority.} \end{aligned} } \]

State feedback assumes that all states are measured. Observer design removes that assumption. A Luenberger observer is

\[ \dot{\hat{\mathbf{x} } }=A\hat{\mathbf{x} }+B\mathbf{u} +L(\mathbf{y}-C\hat{\mathbf{x} }). \]

With estimation error \(\mathbf{e}=\mathbf{x}-\hat{\mathbf{x} }\), and assuming \(D=0\) for simplicity,

\[ \dot{\mathbf{e} }=(A-LC)\mathbf{e}. \]

Thus observer design is dual to state-feedback design: \((A,C)\) observable permits assigning the eigenvalues of \(A-LC\), just as \((A,B)\) controllable permits assigning the eigenvalues of \(A-BK\).

Separation principle. Use

\[ \mathbf{u}=-K\hat{\mathbf{x} }. \]

In coordinates \((\mathbf{x},\mathbf{e})\), the combined dynamics are block upper triangular:

\[ \begin{bmatrix}\dot{\mathbf{x} }\\\dot{\mathbf{e} }\end{bmatrix} = \begin{bmatrix} A-BK & BK \\ 0 & A-LC \end{bmatrix} \begin{bmatrix}\mathbf{x}\\\mathbf{e}\end{bmatrix}. \]

Therefore, the eigenvalues of the combined observer-based controller are

\[ \sigma\!\left( \begin{bmatrix} A-BK & BK \\ 0 & A-LC \end{bmatrix}\right) = \sigma(A-BK)\cup\sigma(A-LC). \]

Controller poles and observer poles can be designed separately for LTI systems under the standard assumptions. This is the conceptual doorway to Kalman filtering and output-feedback optimal control.

The following workflow summarizes how a modern control design naturally expands into the advanced courses covered in this lesson.

flowchart TD
  A["Start with state-space model"] --> B["Check controllability and observability"]
  B --> C["Choose nominal state-feedback design"]
  C --> D["Add quadratic weights for LQR"]
  D --> E["Add estimator if states are not measured"]
  E --> F["Check robustness to uncertainty"]
  F --> G["Add constraints and horizon for MPC"]
  G --> H["Analyze MIMO directions and coupling"]
  H --> I["Validate by simulation and experiments"]
        

The same state-space operations appear across languages, but different ecosystems emphasize different workflows.

• Python: NumPy, SciPy, python-control, cvxpy, OSQP, and Slycot support simulation, Riccati equations, pole placement, frequency response, and convex optimization.
• C++: Eigen, Armadillo, OSQP, qpOASES, CasADi, and ACADO are common for embedded MPC and real-time numerical optimization.
• Java: EJML, Apache Commons Math, and ojAlgo can support educational state-space simulation, matrix computation, and optimization prototypes.
• MATLAB/Simulink: Control System Toolbox, Robust Control Toolbox, Model Predictive Control Toolbox, Simulink, and Simscape are standard for academic and industrial controller prototyping.
• Wolfram Mathematica: symbolic matrix manipulation, exact eigenvalue calculations, and control-system functions are useful for deriving small examples and verifying algebra.

This script compares pole placement, LQR, observer design, a scalar robustness scan, and unconstrained finite-horizon MPC for the same double-integrator plant.

# Chapter30_Lesson5.py
"""
Bridge lab for Modern Control, Chapter 30 Lesson 5.

Topic:
    Connections among pole placement, LQR, robustness checks,
    finite-horizon MPC, MIMO/state-space thinking, and observer design.

Dependencies:
    numpy, scipy
Optional:
    matplotlib for plots if you extend the script.

Run:
    python Chapter30_Lesson5.py
"""
from __future__ import annotations

import numpy as np
from numpy.linalg import inv
from scipy.linalg import solve_continuous_are, eigvals
from scipy.signal import place_poles, cont2discrete


def lqr_gain(A: np.ndarray, B: np.ndarray, Q: np.ndarray, R: np.ndarray):
    """Continuous-time infinite-horizon LQR."""
    P = solve_continuous_are(A, B, Q, R)
    K = inv(R) @ B.T @ P
    closed_loop_poles = eigvals(A - B @ K)
    return K, P, closed_loop_poles


def observer_gain(A: np.ndarray, C: np.ndarray, desired_poles):
    """Luenberger observer gain by dual pole placement."""
    placed = place_poles(A.T, C.T, desired_poles)
    L = placed.gain_matrix.T
    return L, eigvals(A - L @ C)


def pole_placement_gain(A: np.ndarray, B: np.ndarray, desired_poles):
    """State-feedback pole placement."""
    placed = place_poles(A, B, desired_poles)
    K = placed.gain_matrix
    return K, eigvals(A - B @ K)


def robustness_scan(A: np.ndarray, B: np.ndarray, K: np.ndarray, E: np.ndarray, deltas):
    """
    Scan a scalar structured uncertainty A_delta = A + delta E.
    This is not a proof of robust stability; it is a numerical warning test.
    """
    rows = []
    for delta in deltas:
        lam = eigvals(A + delta * E - B @ K)
        rows.append((float(delta), float(np.max(np.real(lam)))))
    return rows


def finite_horizon_lqr_gain(Ad: np.ndarray, Bd: np.ndarray, Q: np.ndarray,
                            R: np.ndarray, Qf: np.ndarray, horizon: int):
    """
    Unconstrained finite-horizon MPC/LQR.
    The first control law is u_0 = -K_0 x_0.
    """
    P = Qf.copy()
    gains = []
    for _ in range(horizon, 0, -1):
        S = R + Bd.T @ P @ Bd
        K = inv(S) @ (Bd.T @ P @ Ad)
        gains.append(K)
        P = Q + Ad.T @ P @ (Ad - Bd @ K)
    gains.reverse()
    return gains[0], gains, P


def simulate_output_feedback(A: np.ndarray, B: np.ndarray, C: np.ndarray,
                             K: np.ndarray, L: np.ndarray,
                             x0: np.ndarray, xhat0: np.ndarray,
                             dt: float = 0.001, tf: float = 4.0):
    """
    Euler simulation of observer-based state feedback:
        u = -K xhat
        xhat_dot = A xhat + B u + L(y - C xhat)
    """
    steps = int(tf / dt)
    x = x0.astype(float).copy()
    xhat = xhat0.astype(float).copy()
    history = []
    for k in range(steps + 1):
        t = k * dt
        y = C @ x
        u = -K @ xhat
        history.append((t, x.copy(), xhat.copy(), float(u[0, 0])))
        xdot = A @ x + B @ u
        xhat_dot = A @ xhat + B @ u + L @ (y - C @ xhat)
        x += dt * xdot
        xhat += dt * xhat_dot
    return history


def main():
    # Double-integrator plant: position and velocity states.
    A = np.array([[0.0, 1.0],
                  [0.0, 0.0]])
    B = np.array([[0.0],
                  [1.0]])
    C = np.array([[1.0, 0.0]])

    # 1) Pole placement: choose closed-loop modes directly.
    K_pp, pp_poles = pole_placement_gain(A, B, [-2.0, -3.0])

    # 2) LQR: choose a quadratic performance index instead of explicit poles.
    Q = np.diag([10.0, 1.0])
    R = np.array([[0.25]])
    K_lqr, P, lqr_poles = lqr_gain(A, B, Q, R)

    # 3) Observer: estimate states from output y = position.
    L, observer_poles = observer_gain(A, C, [-8.0, -9.0])

    # 4) Robustness preview: perturb the spring-like term A[1,0].
    E = np.array([[0.0, 0.0],
                  [1.0, 0.0]])
    deltas = np.linspace(-4.0, 4.0, 17)
    scan = robustness_scan(A, B, K_lqr, E, deltas)

    # 5) MPC preview: discretize and compute first finite-horizon gain.
    Ts = 0.05
    Ad, Bd, Cd, Dd, _ = cont2discrete((A, B, C, np.zeros((1, 1))), Ts)
    K_mpc0, gains, _ = finite_horizon_lqr_gain(
        Ad, Bd, Q=np.diag([10.0, 1.0]), R=np.array([[0.25]]),
        Qf=P, horizon=25
    )

    # 6) Separation principle simulation.
    history = simulate_output_feedback(
        A, B, C, K_lqr, L,
        x0=np.array([[1.0], [0.0]]),
        xhat0=np.array([[0.0], [0.0]]),
        dt=0.002,
        tf=4.0
    )

    print("Pole-placement K:", K_pp)
    print("Pole-placement closed-loop poles:", pp_poles)
    print("\nLQR K:", K_lqr)
    print("LQR Riccati P:\n", P)
    print("LQR closed-loop poles:", lqr_poles)
    print("\nObserver L:\n", L)
    print("Observer poles:", observer_poles)
    print("\nFinite-horizon MPC first gain K0:", K_mpc0)
    print("\nRobustness scan: delta, max real eigenvalue")
    for delta, max_real in scan:
        stable = "stable" if max_real < 0.0 else "unstable"
        print(f"{delta:+.2f}  {max_real:+.4f}  {stable}")

    final_t, final_x, final_xhat, final_u = history[-1]
    print("\nOutput-feedback simulation final state:")
    print("t =", final_t)
    print("x =", final_x.ravel())
    print("xhat =", final_xhat.ravel())
    print("u =", final_u)


if __name__ == "__main__":
    main()

      

This C++ implementation shows the finite-horizon MPC/LQR recursion from scratch for a two-state sampled system, suitable for embedded-control discussion.

// Chapter30_Lesson5.cpp
//
// Bridge lab for Modern Control, Chapter 30 Lesson 5.
// Implements an unconstrained finite-horizon MPC/LQR recursion from scratch
// for a sampled double integrator.
//
// Compile:
//   g++ -std=c++17 -O2 Chapter30_Lesson5.cpp -o Chapter30_Lesson5
//
// Run:
//   ./Chapter30_Lesson5

#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>

struct Vec2 {
    double x1;
    double x2;
};

struct Mat2 {
    double a11, a12, a21, a22;
};

Vec2 add(Vec2 a, Vec2 b) {
    return {a.x1 + b.x1, a.x2 + b.x2};
}

Vec2 sub(Vec2 a, Vec2 b) {
    return {a.x1 - b.x1, a.x2 - b.x2};
}

Vec2 scale(Vec2 a, double s) {
    return {s * a.x1, s * a.x2};
}

Vec2 mat_vec(Mat2 A, Vec2 x) {
    return {A.a11 * x.x1 + A.a12 * x.x2,
            A.a21 * x.x1 + A.a22 * x.x2};
}

Mat2 add(Mat2 A, Mat2 B) {
    return {A.a11 + B.a11, A.a12 + B.a12,
            A.a21 + B.a21, A.a22 + B.a22};
}

Mat2 sub(Mat2 A, Mat2 B) {
    return {A.a11 - B.a11, A.a12 - B.a12,
            A.a21 - B.a21, A.a22 - B.a22};
}

Mat2 mul(Mat2 A, Mat2 B) {
    return {
        A.a11 * B.a11 + A.a12 * B.a21,
        A.a11 * B.a12 + A.a12 * B.a22,
        A.a21 * B.a11 + A.a22 * B.a21,
        A.a21 * B.a12 + A.a22 * B.a22
    };
}

Mat2 transpose(Mat2 A) {
    return {A.a11, A.a21, A.a12, A.a22};
}

double dot(Vec2 a, Vec2 b) {
    return a.x1 * b.x1 + a.x2 * b.x2;
}

Vec2 row_times_mat(Vec2 row, Mat2 A) {
    return {row.x1 * A.a11 + row.x2 * A.a21,
            row.x1 * A.a12 + row.x2 * A.a22};
}

Mat2 outer(Vec2 a, Vec2 b) {
    return {a.x1 * b.x1, a.x1 * b.x2,
            a.x2 * b.x1, a.x2 * b.x2};
}

Vec2 mat_transpose_vec(Mat2 A, Vec2 x) {
    return {A.a11 * x.x1 + A.a21 * x.x2,
            A.a12 * x.x1 + A.a22 * x.x2};
}

struct GainResult {
    Vec2 K0;
    std::vector<Vec2> gains;
};

GainResult finite_horizon_lqr(Mat2 Ad, Vec2 Bd, Mat2 Q, double R,
                              Mat2 Qf, int horizon) {
    Mat2 P = Qf;
    std::vector<Vec2> backward;

    for (int k = horizon - 1; k >= 0; --k) {
        Vec2 PB = mat_vec(P, Bd);
        double S = R + dot(Bd, PB);
        Vec2 PA_col1 = mat_vec(P, {Ad.a11, Ad.a21});
        Vec2 PA_col2 = mat_vec(P, {Ad.a12, Ad.a22});

        Vec2 BtPA = {
            (Bd.x1 * PA_col1.x1 + Bd.x2 * PA_col1.x2) / S,
            (Bd.x1 * PA_col2.x1 + Bd.x2 * PA_col2.x2) / S
        };
        backward.push_back(BtPA);

        Mat2 BK = outer(Bd, BtPA);
        Mat2 Ad_minus_BK = sub(Ad, BK);
        P = add(Q, mul(transpose(Ad), mul(P, Ad_minus_BK)));
    }

    std::vector<Vec2> gains(backward.rbegin(), backward.rend());
    return {gains.front(), gains};
}

int main() {
    const double Ts = 0.05;

    // Sampled double integrator:
    // x_{k+1} = Ad x_k + Bd u_k
    Mat2 Ad{1.0, Ts, 0.0, 1.0};
    Vec2 Bd{0.5 * Ts * Ts, Ts};

    Mat2 Q{10.0, 0.0, 0.0, 1.0};
    double R = 0.25;
    Mat2 Qf{10.0, 0.0, 0.0, 1.0};

    GainResult result = finite_horizon_lqr(Ad, Bd, Q, R, Qf, 25);
    Vec2 K = result.K0;

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Finite-horizon MPC first gain K0 = ["
              << K.x1 << ", " << K.x2 << "]\n";

    Vec2 x{1.0, 0.0};
    double cost = 0.0;

    for (int k = 0; k < 80; ++k) {
        double u = -(K.x1 * x.x1 + K.x2 * x.x2);
        double stage = 10.0 * x.x1 * x.x1 + x.x2 * x.x2 + R * u * u;
        cost += stage;

        Vec2 next = add(mat_vec(Ad, x), scale(Bd, u));
        if (k % 10 == 0) {
            std::cout << "k=" << std::setw(2) << k
                      << "  x=[" << std::setw(10) << x.x1
                      << ", " << std::setw(10) << x.x2
                      << "]  u=" << std::setw(10) << u << "\n";
        }
        x = next;
    }

    std::cout << "Approximate accumulated cost = " << cost << "\n";
    std::cout << "Final state = [" << x.x1 << ", " << x.x2 << "]\n";
    return 0;
}

      

The Java version mirrors the C++ finite-horizon recursion and is useful for students who want to study matrix algorithms without relying on a specialized control toolbox.

// Chapter30_Lesson5.java
//
// Bridge lab for Modern Control, Chapter 30 Lesson 5.
// Implements an unconstrained finite-horizon MPC/LQR recursion from scratch
// for a sampled double integrator.
//
// Compile:
//   javac Chapter30_Lesson5.java
//
// Run:
//   java Chapter30_Lesson5

import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

public class Chapter30_Lesson5 {
    static class Vec2 {
        double x1, x2;
        Vec2(double x1, double x2) {
            this.x1 = x1;
            this.x2 = x2;
        }
    }

    static class Mat2 {
        double a11, a12, a21, a22;
        Mat2(double a11, double a12, double a21, double a22) {
            this.a11 = a11;
            this.a12 = a12;
            this.a21 = a21;
            this.a22 = a22;
        }
    }

    static class GainResult {
        Vec2 K0;
        List<Vec2> gains;
        GainResult(Vec2 K0, List<Vec2> gains) {
            this.K0 = K0;
            this.gains = gains;
        }
    }

    static Vec2 add(Vec2 a, Vec2 b) {
        return new Vec2(a.x1 + b.x1, a.x2 + b.x2);
    }

    static Vec2 scale(Vec2 a, double s) {
        return new Vec2(s * a.x1, s * a.x2);
    }

    static Vec2 matVec(Mat2 A, Vec2 x) {
        return new Vec2(
            A.a11 * x.x1 + A.a12 * x.x2,
            A.a21 * x.x1 + A.a22 * x.x2
        );
    }

    static Mat2 add(Mat2 A, Mat2 B) {
        return new Mat2(
            A.a11 + B.a11, A.a12 + B.a12,
            A.a21 + B.a21, A.a22 + B.a22
        );
    }

    static Mat2 sub(Mat2 A, Mat2 B) {
        return new Mat2(
            A.a11 - B.a11, A.a12 - B.a12,
            A.a21 - B.a21, A.a22 - B.a22
        );
    }

    static Mat2 mul(Mat2 A, Mat2 B) {
        return new Mat2(
            A.a11 * B.a11 + A.a12 * B.a21,
            A.a11 * B.a12 + A.a12 * B.a22,
            A.a21 * B.a11 + A.a22 * B.a21,
            A.a21 * B.a12 + A.a22 * B.a22
        );
    }

    static Mat2 transpose(Mat2 A) {
        return new Mat2(A.a11, A.a21, A.a12, A.a22);
    }

    static double dot(Vec2 a, Vec2 b) {
        return a.x1 * b.x1 + a.x2 * b.x2;
    }

    static Mat2 outer(Vec2 a, Vec2 b) {
        return new Mat2(
            a.x1 * b.x1, a.x1 * b.x2,
            a.x2 * b.x1, a.x2 * b.x2
        );
    }

    static GainResult finiteHorizonLqr(Mat2 Ad, Vec2 Bd, Mat2 Q, double R,
                                       Mat2 Qf, int horizon) {
        Mat2 P = Qf;
        List<Vec2> backward = new ArrayList<>();

        for (int k = horizon - 1; k >= 0; --k) {
            Vec2 PB = matVec(P, Bd);
            double S = R + dot(Bd, PB);

            Vec2 PAcol1 = matVec(P, new Vec2(Ad.a11, Ad.a21));
            Vec2 PAcol2 = matVec(P, new Vec2(Ad.a12, Ad.a22));

            Vec2 K = new Vec2(
                (Bd.x1 * PAcol1.x1 + Bd.x2 * PAcol1.x2) / S,
                (Bd.x1 * PAcol2.x1 + Bd.x2 * PAcol2.x2) / S
            );
            backward.add(K);

            Mat2 BK = outer(Bd, K);
            Mat2 AdMinusBK = sub(Ad, BK);
            P = add(Q, mul(transpose(Ad), mul(P, AdMinusBK)));
        }

        Collections.reverse(backward);
        return new GainResult(backward.get(0), backward);
    }

    public static void main(String[] args) {
        double Ts = 0.05;

        Mat2 Ad = new Mat2(1.0, Ts, 0.0, 1.0);
        Vec2 Bd = new Vec2(0.5 * Ts * Ts, Ts);

        Mat2 Q = new Mat2(10.0, 0.0, 0.0, 1.0);
        double R = 0.25;
        Mat2 Qf = new Mat2(10.0, 0.0, 0.0, 1.0);

        GainResult result = finiteHorizonLqr(Ad, Bd, Q, R, Qf, 25);
        Vec2 K = result.K0;

        System.out.printf("Finite-horizon MPC first gain K0 = [%.6f, %.6f]%n",
                          K.x1, K.x2);

        Vec2 x = new Vec2(1.0, 0.0);
        double cost = 0.0;

        for (int k = 0; k < 80; ++k) {
            double u = -(K.x1 * x.x1 + K.x2 * x.x2);
            double stage = 10.0 * x.x1 * x.x1 + x.x2 * x.x2 + R * u * u;
            cost += stage;

            if (k % 10 == 0) {
                System.out.printf("k=%2d  x=[%10.6f, %10.6f]  u=%10.6f%n",
                                  k, x.x1, x.x2, u);
            }

            x = add(matVec(Ad, x), scale(Bd, u));
        }

        System.out.printf("Approximate accumulated cost = %.6f%n", cost);
        System.out.printf("Final state = [%.6f, %.6f]%n", x.x1, x.x2);
    }
}

      

The MATLAB file uses standard design commands where available and also includes a local finite-horizon Riccati recursion. If Simulink is installed, it creates a minimal state-space model skeleton.

% Chapter30_Lesson5.m
% Bridge lab for Modern Control, Chapter 30 Lesson 5.
%
% Requires for full functionality:
%   Control System Toolbox
% Optional:
%   Simulink, Robust Control Toolbox, Model Predictive Control Toolbox
%
% Run:
%   Chapter30_Lesson5

clear; clc;

A = [0 1; 0 0];
B = [0; 1];
C = [1 0];
D = 0;

Q = diag([10 1]);
R = 0.25;

fprintf('--- Pole placement ---\n');
Kpp = place(A, B, [-2 -3]);
disp('Kpp ='); disp(Kpp);
disp('eig(A-B*Kpp) ='); disp(eig(A-B*Kpp).');

fprintf('\n--- LQR ---\n');
[Klqr, P, polesLqr] = lqr(A, B, Q, R);
disp('Klqr ='); disp(Klqr);
disp('Riccati matrix P ='); disp(P);
disp('LQR closed-loop poles ='); disp(polesLqr.');

fprintf('\n--- Observer design by duality ---\n');
L = place(A', C', [-8 -9])';
disp('L ='); disp(L);
disp('eig(A-L*C) ='); disp(eig(A-L*C).');

fprintf('\n--- Separation principle check ---\n');
Asep = [A-B*Klqr, B*Klqr; zeros(2), A-L*C];
disp('eig(separation matrix) ='); disp(eig(Asep).');

fprintf('\n--- Robustness preview: scalar structured uncertainty ---\n');
E = [0 0; 1 0];
for delta = -4:0.5:4
    lam = eig(A + delta*E - B*Klqr);
    maxReal = max(real(lam));
    if maxReal < 0
        status = 'stable';
    else
        status = 'unstable';
    end
    fprintf('delta=%+5.2f  max_real=%+8.4f  %s\n', delta, maxReal, status);
end

fprintf('\n--- Finite-horizon unconstrained MPC/LQR ---\n');
Ts = 0.05;
sysd = c2d(ss(A,B,C,D), Ts);
Ad = sysd.A;
Bd = sysd.B;

Qd = Q;
Rd = R;
Qf = P;
N = 25;

[Pseq, Kseq] = finite_horizon_lqr_local(Ad, Bd, Qd, Rd, Qf, N);
K0 = Kseq(:,:,1);
disp('First MPC gain K0 ='); disp(K0);

x = [1; 0];
cost = 0;
for k = 1:80
    u = -K0*x;
    cost = cost + x'*Qd*x + u'*Rd*u;
    x = Ad*x + Bd*u;
end
disp('Final state after repeated first-gain simulation ='); disp(x.');
fprintf('Approximate accumulated cost = %.6f\n', cost);

fprintf('\n--- Optional Simulink skeleton ---\n');
if license('test','Simulink')
    model = 'Chapter30_Lesson5_Simulink_Observer_LQR';
    if bdIsLoaded(model), close_system(model, 0); end
    new_system(model);
    open_system(model);
    add_block('simulink/Sources/Step', [model '/Initial command']);
    add_block('simulink/Continuous/State-Space', [model '/Plant']);
    add_block('simulink/Sinks/Scope', [model '/Scope']);
    set_param([model '/Plant'], 'A', 'A-B*Klqr', 'B', 'B', 'C', 'C', 'D', 'D');
    add_line(model, 'Initial command/1', 'Plant/1');
    add_line(model, 'Plant/1', 'Scope/1');
    save_system(model);
    fprintf('Created Simulink model: %s.slx\n', model);
else
    fprintf('Simulink license not detected; skipped model generation.\n');
end

function [Pseq, Kseq] = finite_horizon_lqr_local(A, B, Q, R, Qf, N)
    nx = size(A,1);
    nu = size(B,2);
    Pseq = zeros(nx,nx,N+1);
    Kseq = zeros(nu,nx,N);
    Pseq(:,:,N+1) = Qf;
    for k = N:-1:1
        Pnext = Pseq(:,:,k+1);
        S = R + B'*Pnext*B;
        K = S\(B'*Pnext*A);
        Kseq(:,:,k) = K;
        Pseq(:,:,k) = Q + A'*Pnext*(A-B*K);
    end
end

      

The Mathematica notebook emphasizes exact symbolic setup for the Riccati equation and direct verification of observer and MPC calculations.

Notebook[{
Cell["Chapter30_Lesson5.nb", "Title"],
Cell["Bridge lab for Modern Control: LQR, observer design, robustness preview, and finite-horizon MPC.", "Text"],
Cell[BoxData[
RowBox[{
"ClearAll", "[", "\"Global`*\"", "]"}]], "Input"],
Cell[BoxData[
RowBox[{
RowBox[{"A", "=", "{ {0,1},{0,0} }"}], ";",
RowBox[{"B", "=", "{ {0},{1} }"}], ";",
RowBox[{"Cmat", "=", "{ {1,0} }"}], ";",
RowBox[{"Q", "=", "DiagonalMatrix[{10,1}]"}], ";",
RowBox[{"R", "=", "{ {1/4} }"}], ";"}]], "Input"],
Cell["Continuous-time LQR from the algebraic Riccati equation.", "Section"],
Cell[BoxData[
RowBox[{
RowBox[{"P", "=", "{ {p11,p12},{p12,p22} }"}], ";",
RowBox[{"care", "=", 
RowBox[{
RowBox[{"Transpose[A].P"}], "+", RowBox[{"P.A"}], "-", 
RowBox[{"P.B.Inverse[R].Transpose[B].P"}], "+", "Q"}]}], ";",
RowBox[{"sols", "=", RowBox[{"Solve", "[", 
RowBox[{RowBox[{"Thread", "[", RowBox[{"Flatten[care]", "==", "0"}], "]"}], ",", 
"{p11,p12,p22}", ",", "Reals"}], "]"}]}], ";",
RowBox[{"Psol", "=", RowBox[{"N", "[", RowBox[{"P", "/.", RowBox[{"sols", "[[4]]"}]}], "]"}]}], ";",
RowBox[{"Klqr", "=", RowBox[{"N", "[", RowBox[{"Inverse[R].Transpose[B].Psol"}], "]"}]}], ";",
RowBox[{"Eigenvalues", "[", RowBox[{"A", "-", RowBox[{"B.Klqr"}]}], "]"}]}]], "Input"],
Cell["Observer gain by direct characteristic-polynomial matching for the double integrator.", "Section"],
Cell[BoxData[
RowBox[{
RowBox[{"L", "=", "{ {17},{72} }"}], ";",
RowBox[{"Eigenvalues", "[", RowBox[{"A", "-", RowBox[{"L.Cmat"}]}], "]"}]}]], "Input"],
Cell["Finite-horizon unconstrained MPC recursion for a sampled double integrator.", "Section"],
Cell[BoxData[
RowBox[{
RowBox[{"Ts", "=", "0.05"}], ";",
RowBox[{"Ad", "=", "{ {1,Ts},{0,1} }"}], ";",
RowBox[{"Bd", "=", "{ {Ts^2/2},{Ts} }"}], ";",
RowBox[{"Qd", "=", "Q"}], ";",
RowBox[{"Rd", "=", "R"}], ";",
RowBox[{"Nhor", "=", "25"}], ";",
RowBox[{"Plist", "=", "{Psol}"}], ";",
RowBox[{"Klist", "=", "{}"}], ";",
RowBox[{"Do", "[", 
RowBox[{
RowBox[{
RowBox[{"Pnext", "=", RowBox[{"First", "[", "Plist", "]"}]}], ";",
RowBox[{"S", "=", RowBox[{"Rd", "+", RowBox[{"Transpose[Bd].Pnext.Bd"}]}]}], ";",
RowBox[{"K", "=", RowBox[{"Inverse[S].Transpose[Bd].Pnext.Ad"}]}], ";",
RowBox[{"PrependTo", "[", RowBox[{"Klist", ",", "K"}], "]"}], ";",
RowBox[{"PrependTo", "[", RowBox[{"Plist", ",", 
RowBox[{"Qd", "+", RowBox[{"Transpose[Ad].Pnext.", RowBox[{"(", RowBox[{"Ad", "-", RowBox[{"Bd.K"}]}], ")"}]}]}]}], "]"}]}],
",", "{Nhor}"}], "]"}], ";",
RowBox[{"First", "[", "Klist", "]"}]}]], "Input"],
Cell["Robustness scan for A_delta = A + delta E.", "Section"],
Cell[BoxData[
RowBox[{
RowBox[{"Eunc", "=", "{ {0,0},{1,0} }"}], ";",
RowBox[{"Table", "[", 
RowBox[{
RowBox[{"{", RowBox[{"delta", ",", 
RowBox[{"Max", "[", RowBox[{"Re", "[", RowBox[{"Eigenvalues", "[", 
RowBox[{"A", "+", RowBox[{"delta Eunc"}], "-", RowBox[{"B.Klqr"}]}], "]"}], "]"}], "]"}]}], "}"}],
",", RowBox[{"{", RowBox[{"delta", ",", "-4", ",", "4", ",", "0.5"}], "}"}]}], "]"}]}]], "Input"]
}]

      

Problem 1 (LQR as a stabilizing state-feedback law): Suppose \(P=P^\top\succeq 0\) solves \(A^\top P+PA-PBR^{-1}B^\top P+Q=0\) with \(R\succ 0\). Show that \(K=R^{-1}B^\top P\) gives \(\dot{V}\leq 0\) for \(V=\mathbf{x}^\top P\mathbf{x}\).

Solution: Substituting \(\mathbf{u}=-K\mathbf{x}\) gives \(\dot{\mathbf{x} }=(A-BK)\mathbf{x}\). Therefore,

\[ \dot{V} =\mathbf{x}^\top((A-BK)^\top P+P(A-BK))\mathbf{x} =-\mathbf{x}^\top(Q+K^\top RK)\mathbf{x}\leq 0. \]

Since \(Q\succeq 0\) and \(R\succ 0\), the expression is nonpositive. With the usual stabilizability and detectability assumptions, the stabilizing Riccati solution makes \(A-BK\) Hurwitz.

Problem 2 (Separation principle): For \(\mathbf{u}=-K\hat{\mathbf{x} }\) and observer \(\dot{\hat{\mathbf{x} } }=A\hat{\mathbf{x} }+B\mathbf{u} +L(\mathbf{y}-C\hat{\mathbf{x} })\), prove that controller and observer eigenvalues appear independently in the combined system.

Solution: Let \(\mathbf{e}=\mathbf{x}-\hat{\mathbf{x} }\). Then \(\dot{\mathbf{e} }=(A-LC)\mathbf{e}\). In \((\mathbf{x},\mathbf{e})\) coordinates,

\[ \begin{bmatrix}\dot{\mathbf{x} }\\\dot{\mathbf{e} }\end{bmatrix} = \begin{bmatrix}A-BK & BK\\0&A-LC\end{bmatrix} \begin{bmatrix}\mathbf{x}\\\mathbf{e}\end{bmatrix}. \]

The matrix is block upper triangular, so its characteristic polynomial is the product of the characteristic polynomials of \(A-BK\) and \(A-LC\).

Problem 3 (MPC without constraints): Show that unconstrained finite-horizon MPC is equivalent to a finite-horizon LQR recursion.

Solution: Assume the value function at stage \(k+1\) is \(V_{k+1}(\mathbf{x})=\mathbf{x}^\top P_{k+1}\mathbf{x}\). The one-step dynamic-programming minimization is

\[ \min_{\mathbf{u}_k} \left[ \mathbf{x}_k^\top Q\mathbf{x}_k+\mathbf{u}_k^\top R\mathbf{u}_k+ (A_d\mathbf{x}_k+B_d\mathbf{u}_k)^\top P_{k+1} (A_d\mathbf{x}_k+B_d\mathbf{u}_k) \right]. \]

Differentiating with respect to \(\mathbf{u}_k\) gives the linear optimal law \(\mathbf{u}_k=-K_k\mathbf{x}_k\), where

\[ K_k=(R+B_d^\top P_{k+1}B_d)^{-1}B_d^\top P_{k+1}A_d. \]

Substitution gives the Riccati recursion for \(P_k\).

Problem 4 (Robust stability by a common Lyapunov matrix): Consider a family \(\mathcal{A}=\{A_1,A_2,\ldots,A_m\}\). State a sufficient condition that guarantees all systems \(\dot{\mathbf{x} }=A_i\mathbf{x}\) are stable.

Solution: If there exists \(P=P^\top\succ 0\) such that

\[ A_i^\top P+PA_i\prec 0,\qquad i=1,2,\ldots,m, \]

then \(V=\mathbf{x}^\top P\mathbf{x}\) decreases for every member of the family. Hence all systems in the family are asymptotically stable under the same quadratic certificate.

Problem 5 (MIMO directional gain): Let \(\mathbf{y}=G(j\omega)\mathbf{u}\). Show why the largest singular value is the maximum amplification over all unit-norm input directions.

Solution: For \(\|\mathbf{u}\|_2=1\),

\[ \|\mathbf{y}\|_2^2 =\mathbf{u}^*G(j\omega)^*G(j\omega)\mathbf{u}. \]

By the Rayleigh-Ritz theorem, the maximum over unit vectors is \(\lambda_{\max}(G^*G)\). Taking the square root gives \(\bar{\sigma}(G(j\omega))\).

Optimal control chooses feedback through a cost functional. Robust control asks whether the design survives uncertainty. MPC repeatedly solves a finite-horizon constrained state-space optimization problem. Multivariable control studies directions, coupling, and interaction. Observer design reconstructs missing states and enables output feedback. All of these courses are natural continuations of the state-space concepts developed throughout Modern Control.