In the previous two lessons, quadratic forms and signal-energy ideas were introduced as ways of assigning a scalar measure to vectors and signals. Here we specialize those ideas to a continuous-time state-space system with state \( x(t)\in\mathbb{R}^n \) and input \( u(t)\in\mathbb{R}^m \). A weighted performance descriptor over a finite horizon \( [0,T] \) is

\[ J_T(x,u)=\int_0^T \left(x(t)^\top Qx(t)+u(t)^\top Ru(t)\right)dt. \]

The matrix \( Q \) assigns relative importance to state deviations, while \( R \) assigns relative importance to control usage. At this stage, we do not yet solve an optimal-control problem. Instead, we treat \( Q \) and \( R \) as precise mathematical descriptions of what a designer means by a “large state error” and a “large input effort.”

flowchart TD
  A["Engineering priorities"] --> B["State limits and state importance"]
  A --> C["Actuator limits and input cost"]
  B --> D["Choose Q"]
  C --> E["Choose R"]
  D --> F["Weighted state term: xT Q x"]
  E --> G["Weighted input term: uT R u"]
  F --> H["Performance descriptor J"]
  G --> H
  H --> I["Compare trajectories and feedback designs"]
        

For real-valued quadratic measures, the usual assumptions are \( Q=Q^\top\succeq 0 \) and \( R=R^\top\succ 0 \). The positive semidefinite requirement on \( Q \) permits some state directions to be unpenalized. The positive definite requirement on \( R \) means every nonzero input has a strictly positive instantaneous cost.

\[ Q\in\mathbb{S}_+^n,\qquad R\in\mathbb{S}_{++}^m,\qquad x^\top Qx\ge 0,\qquad u^\top Ru>0\;\text{for}\;u\ne 0. \]

The expression above uses the cone notation \( \mathbb{S}_+^n \) for symmetric positive semidefinite matrices and \( \mathbb{S}_{++}^m \) for symmetric positive definite matrices. If a designer writes a nonsymmetric matrix in a quadratic form, only its symmetric part matters:

\[ x^\top Qx=x^\top\left(\frac{Q+Q^\top}{2}\right)x. \]

Proof. Decompose \( Q \) into its symmetric and skew-symmetric parts:

\[ Q=\frac{Q+Q^\top}{2}+\frac{Q-Q^\top}{2}=Q_s+Q_a, \qquad Q_a^\top=-Q_a. \]

For any real vector \( x \), the scalar \( x^\top Q_ax \) equals its own transpose:

\[ x^\top Q_ax=(x^\top Q_ax)^\top=x^\top Q_a^\top x=-x^\top Q_ax. \]

Hence \( x^\top Q_ax=0 \), proving that the skew-symmetric part has no effect on the quadratic performance measure.

If \( Q\succ 0 \), the set of states with equal weighted magnitude is an ellipsoid:

\[ \mathcal{E}_Q(c)=\left\{x\in\mathbb{R}^n\;:\;x^\top Qx=c\right\}. \]

Let the spectral decomposition of \( Q \) be \( Q=V\Lambda V^\top \), with orthonormal eigenvectors in \( V \) and nonnegative eigenvalues in \( \Lambda \). In modal coordinates \( z=V^\top x \), the quadratic form becomes

\[ x^\top Qx=z^\top\Lambda z=\sum_{i=1}^n \lambda_i z_i^2. \]

Large eigenvalues compress the ellipsoid along the corresponding directions, meaning the performance descriptor regards those directions as expensive. A diagonal \( Q \) penalizes individual coordinates independently, while a nondiagonal \( Q \) introduces cross-weighting between state coordinates.

The same interpretation applies to \( R \) in the input space. For multi-input systems, off-diagonal terms in \( R \) can represent coupled actuator usage, shared power supplies, or a preference against simultaneous use of certain channels.

The state vector may contain variables with different physical units: position, velocity, current, angle, angular velocity, pressure, or temperature. A good \( Q \) should therefore be built from physical performance variables rather than arbitrary coordinate magnitudes. Suppose the variables of interest are \( z=C_zx \), where \( z \) collects physically meaningful errors. If \( W_z\succeq 0 \) weights those errors, then

\[ z^\top W_zz=(C_zx)^\top W_z(C_zx)=x^\top(C_z^\top W_zC_z)x. \]

Therefore an induced state weighting matrix is

\[ Q=C_z^\top W_zC_z. \]

This construction is useful when the raw state is not itself the performance output. For example, in a mechanical system, the state might contain both positions and velocities, while the designer might care mainly about selected displacements and a limited set of velocities.

A common first-pass normalization is Bryson-style weighting. If \( x_i^{\max} \) is the largest acceptable magnitude of state component \( x_i \), then the diagonal choice

\[ Q=\operatorname{diag}\left(\frac{1}{(x_1^{\max})^2}, \frac{1}{(x_2^{\max})^2},\dots, \frac{1}{(x_n^{\max})^2}\right) \]

makes a state component near its acceptable limit contribute approximately one unit to \( x^\top Qx \).

The matrix \( R \) describes how expensive control signals are. If all actuators are independent and the acceptable input magnitudes are \( u_j^{\max} \), a normalized diagonal choice is

\[ R=\operatorname{diag}\left(\frac{1}{(u_1^{\max})^2}, \frac{1}{(u_2^{\max})^2},\dots, \frac{1}{(u_m^{\max})^2}\right). \]

If actuator channels are coupled, use a full symmetric positive definite matrix. For example, with two inputs,

\[ u^\top Ru=r_1u_1^2+2r_{12}u_1u_2+r_2u_2^2, \qquad R=\begin{bmatrix}r_1&r_{12}\\r_{12}&r_2\end{bmatrix}\succ 0. \]

The cross-term coefficient \( r_{12} \) changes the cost of using the two inputs together. Positive definiteness is essential: otherwise the performance descriptor could assign zero or negative cost to a nonzero actuator command, which is physically inappropriate for an effort penalty.

State coordinates in modern control are not unique. If \( x=Tz \), with nonsingular \( T \), the same physical quadratic state measure must satisfy

\[ x^\top Q_xx=z^\top Q_zz. \]

Substituting \( x=Tz \) gives

\[ (Tz)^\top Q_x(Tz)=z^\top(T^\top Q_xT)z. \]

Therefore the coordinate-transformed weighting matrix is

\[ Q_z=T^\top Q_xT. \]

This formula is important because a diagonal matrix in one coordinate system does not generally remain diagonal after a change of coordinates. Arbitrary diagonal choices can therefore unintentionally weight coordinate artifacts rather than physical performance variables.

If the input is transformed by \( u=Sv \), then the corresponding input weighting becomes

\[ R_v=S^\top R_uS. \]

Consider two feedback gains \( K_1 \) and \( K_2 \) applied to the same plant \( \dot{x}=Ax+Bu \) through \( u=-Kx \). For the same initial condition, the closed-loop trajectories satisfy

\[ \dot{x}_i=(A-BK_i)x_i, \qquad u_i=-K_ix_i, \qquad i=1,2. \]

The performance descriptor ranks the two designs through

\[ J_T(K_i)=\int_0^T\left(x_i(t)^\top Qx_i(t)+u_i(t)^\top Ru_i(t)\right)dt. \]

It is crucial that \( Q \) and \( R \) remain fixed during the comparison. Changing the weights changes the meaning of performance itself, so it cannot be used to make a fair comparison between designs.

flowchart TD
  A["Start with physical limits"] --> B["Normalize states and inputs"]
  B --> C["Select initial Q and R"]
  C --> D["Simulate candidate feedback designs"]
  D --> E["Compute same J for every design"]
  E --> F["Inspect state peaks and input peaks"]
  F --> G["Weights reflect priorities?"]
  G -->|"yes"| H["Use descriptor for comparison"]
  G -->|"no"| I["Revise Q or R from engineering requirements"]
  I --> C
        

The instantaneous function \( \ell(x,u)=x^\top Qx+u^\top Ru \) is convex in \( (x,u) \) when \( Q\succeq 0 \) and \( R\succeq 0 \). If \( R\succ 0 \), it is strictly convex in \( u \). In block form,

\[ \ell(x,u)=\begin{bmatrix}x\\u\end{bmatrix}^\top \begin{bmatrix}Q&0\\0&R\end{bmatrix} \begin{bmatrix}x\\u\end{bmatrix}. \]

Since a block-diagonal matrix with positive semidefinite diagonal blocks is positive semidefinite, \( \ell \) is a convex quadratic function. Moreover, if \( \lambda_{\min}(R) \) is the smallest eigenvalue of \( R \), then

\[ u^\top Ru\ge \lambda_{\min}(R)\|u\|_2^2. \]

This lower bound shows why a positive definite \( R \) prevents arbitrarily large input actions from being cost-free.

This script checks definiteness, constructs Bryson-style weights, simulates a closed-loop trajectory, evaluates the finite-horizon weighted cost, and verifies coordinate-transformation consistency.

Chapter28_Lesson3.py

"""
Chapter28_Lesson3.py
Modern Control — Chapter 28, Lesson 3
State and control weighting matrices Q and R as performance descriptors.

This script demonstrates:
1. Symmetry and definiteness checks for Q and R.
2. Bryson-style diagonal weight selection.
3. Finite-horizon weighted performance integral for a simulated state/input record.
4. Coordinate-transformation consistency of a quadratic state cost.
"""

from __future__ import annotations

import numpy as np
from numpy.linalg import eigvalsh
from scipy.integrate import solve_ivp


def symmetrize(M: np.ndarray) -> np.ndarray:
    """Return the symmetric part of a square matrix."""
    return 0.5 * (M + M.T)


def is_psd(M: np.ndarray, tol: float = 1e-10) -> bool:
    """Check positive semidefiniteness through symmetric eigenvalues."""
    Ms = symmetrize(M)
    return bool(np.min(eigvalsh(Ms)) >= -tol)


def is_pd(M: np.ndarray, tol: float = 1e-10) -> bool:
    """Check positive definiteness through symmetric eigenvalues."""
    Ms = symmetrize(M)
    return bool(np.min(eigvalsh(Ms)) > tol)


def bryson_weights(max_abs_values: np.ndarray) -> np.ndarray:
    """
    Bryson-style diagonal weighting:
        weight_i = 1 / allowed_i^2.
    The vector contains acceptable maximum magnitudes of variables.
    """
    max_abs_values = np.asarray(max_abs_values, dtype=float)
    if np.any(max_abs_values <= 0.0):
        raise ValueError("All allowed magnitudes must be positive.")
    return np.diag(1.0 / (max_abs_values**2))


def simulate_closed_loop(A: np.ndarray, B: np.ndarray, K: np.ndarray, x0: np.ndarray, tf: float = 8.0):
    """Simulate xdot = (A - B K) x and u = -K x."""
    Acl = A - B @ K

    def rhs(_t, x):
        return Acl @ x

    sol = solve_ivp(rhs, (0.0, tf), x0, dense_output=False, max_step=0.01, rtol=1e-8, atol=1e-10)
    X = sol.y.T
    U = -(K @ X.T).T
    return sol.t, X, U


def weighted_cost(t: np.ndarray, X: np.ndarray, U: np.ndarray, Q: np.ndarray, R: np.ndarray) -> float:
    """Approximate integral of x'Qx + u'Ru using the trapezoidal rule."""
    state_terms = np.einsum("bi,ij,bj->b", X, Q, X)
    input_terms = np.einsum("bi,ij,bj->b", U, R, U)
    return float(np.trapz(state_terms + input_terms, t))


def main() -> None:
    # A stable feedback example for a second-order plant.
    A = np.array([[0.0, 1.0], [-2.0, -0.4]])
    B = np.array([[0.0], [1.0]])
    K = np.array([[3.0, 2.2]])
    x0 = np.array([1.0, 0.0])

    # Performance-descriptor weights from engineering tolerances.
    # State 1 allowed magnitude: 1.0; state 2 allowed magnitude: 2.0.
    Q = bryson_weights(np.array([1.0, 2.0]))
    # Input allowed magnitude: 0.5.
    R = bryson_weights(np.array([0.5]))

    print("Q =\n", Q)
    print("R =\n", R)
    print("Q is PSD:", is_psd(Q))
    print("R is PD:", is_pd(R))

    t, X, U = simulate_closed_loop(A, B, K, x0)
    J = weighted_cost(t, X, U, Q, R)
    print(f"Finite-horizon weighted cost J_T = {J:.6f}")

    # Coordinate transformation: x = T z.
    T = np.array([[2.0, 0.5], [0.0, 1.5]])
    Z = np.linalg.solve(T, X.T).T
    Qz = T.T @ Q @ T
    J_state_x = float(np.trapz(np.einsum("bi,ij,bj->b", X, Q, X), t))
    J_state_z = float(np.trapz(np.einsum("bi,ij,bj->b", Z, Qz, Z), t))
    print(f"State cost in x-coordinates = {J_state_x:.6f}")
    print(f"State cost in z-coordinates = {J_state_z:.6f}")
    print("Coordinate consistency error:", abs(J_state_x - J_state_z))


if __name__ == "__main__":
    main()

The C++ example uses only the standard library and implements RK4 integration from scratch for a two-state closed-loop system.

Chapter28_Lesson3.cpp

/*
Chapter28_Lesson3.cpp
Modern Control — Chapter 28, Lesson 3
State and control weighting matrices Q and R as performance descriptors.

This standalone C++ example computes a finite-horizon weighted cost
for xdot = (A - B K)x, u = -Kx using RK4 integration.
*/

#include <array>
#include <cmath>
#include <iostream>
#include <stdexcept>

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

Vec2 mat_vec(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]};
}

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

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

double quad2(const Vec2& x, const Mat2& Q) {
    return x[0] * (Q[0][0] * x[0] + Q[0][1] * x[1]) +
           x[1] * (Q[1][0] * x[0] + Q[1][1] * x[1]);
}

double input_value(const std::array<double, 2>& K, const Vec2& x) {
    return -(K[0] * x[0] + K[1] * x[1]);
}

Vec2 rhs(const Mat2& Acl, const Vec2& x) {
    return mat_vec(Acl, x);
}

Vec2 rk4_step(const Mat2& Acl, const Vec2& x, double h) {
    Vec2 k1 = rhs(Acl, x);
    Vec2 k2 = rhs(Acl, add(x, scale(0.5 * h, k1)));
    Vec2 k3 = rhs(Acl, add(x, scale(0.5 * h, k2)));
    Vec2 k4 = rhs(Acl, add(x, scale(h, k3)));
    return add(x, scale(h / 6.0, add(add(k1, scale(2.0, k2)), add(scale(2.0, k3), k4))));
}

bool is_pd_1x1(double r) {
    return r > 0.0;
}

bool is_psd_2x2(const Mat2& Q) {
    // Sylvester condition for 2x2 symmetric PSD:
    // q11 >= 0, q22 >= 0, det(Q) >= 0.
    const double det = Q[0][0] * Q[1][1] - Q[0][1] * Q[1][0];
    return Q[0][0] >= -1e-12 && Q[1][1] >= -1e-12 && det >= -1e-12;
}

int main() {
    // Plant and feedback: xdot = (A - B K)x.
    Mat2 A = { { {0.0, 1.0}, {-2.0, -0.4} } };
    std::array<double, 2> B = {0.0, 1.0};
    std::array<double, 2> K = {3.0, 2.2};

    Mat2 Acl = A;
    Acl[0][0] -= B[0] * K[0];
    Acl[0][1] -= B[0] * K[1];
    Acl[1][0] -= B[1] * K[0];
    Acl[1][1] -= B[1] * K[1];

    // Bryson-style tolerances: x1_max=1, x2_max=2, u_max=0.5.
    Mat2 Q = { { {1.0, 0.0}, {0.0, 1.0 / 4.0} } };
    double R = 1.0 / (0.5 * 0.5);

    if (!is_psd_2x2(Q)) {
        throw std::runtime_error("Q must be positive semidefinite.");
    }
    if (!is_pd_1x1(R)) {
        throw std::runtime_error("R must be positive definite.");
    }

    Vec2 x = {1.0, 0.0};
    double h = 0.001;
    double tf = 8.0;
    int steps = static_cast<int>(tf / h);
    double J = 0.0;

    auto integrand = [&](const Vec2& xv) {
        double u = input_value(K, xv);
        return quad2(xv, Q) + R * u * u;
    };

    double f_old = integrand(x);
    for (int k = 0; k < steps; ++k) {
        Vec2 x_next = rk4_step(Acl, x, h);
        double f_new = integrand(x_next);
        J += 0.5 * h * (f_old + f_new);
        x = x_next;
        f_old = f_new;
    }

    std::cout << "Q = [[1, 0], [0, 0.25]]\n";
    std::cout << "R = " << R << "\n";
    std::cout << "Finite-horizon weighted cost J_T = " << J << "\n";
    return 0;
}

The Java implementation mirrors the C++ version and keeps all linear algebra explicit so students can see how the quadratic cost is computed.

Chapter28_Lesson3.java

/*
Chapter28_Lesson3.java
Modern Control — Chapter 28, Lesson 3
State and control weighting matrices Q and R as performance descriptors.

This standalone Java example computes a finite-horizon weighted cost
for xdot = (A - B K)x, u = -Kx using RK4 integration.
*/

public class Chapter28_Lesson3 {
    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]
        };
    }

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

    static double[] scale(double c, double[] x) {
        return new double[] {c * x[0], c * x[1]};
    }

    static double[] rhs(double[][] Acl, double[] x) {
        return matVec(Acl, x);
    }

    static double[] rk4Step(double[][] Acl, double[] x, double h) {
        double[] k1 = rhs(Acl, x);
        double[] k2 = rhs(Acl, add(x, scale(0.5 * h, k1)));
        double[] k3 = rhs(Acl, add(x, scale(0.5 * h, k2)));
        double[] k4 = rhs(Acl, add(x, scale(h, k3)));
        return add(x, scale(h / 6.0, add(add(k1, scale(2.0, k2)), add(scale(2.0, k3), k4))));
    }

    static double quad(double[] x, double[][] Q) {
        return x[0] * (Q[0][0] * x[0] + Q[0][1] * x[1])
             + x[1] * (Q[1][0] * x[0] + Q[1][1] * x[1]);
    }

    static double input(double[] K, double[] x) {
        return -(K[0] * x[0] + K[1] * x[1]);
    }

    static boolean isPsd2x2(double[][] Q) {
        double det = Q[0][0] * Q[1][1] - Q[0][1] * Q[1][0];
        return Q[0][0] >= -1e-12 && Q[1][1] >= -1e-12 && det >= -1e-12;
    }

    public static void main(String[] args) {
        double[][] A = { {0.0, 1.0}, {-2.0, -0.4} };
        double[] B = {0.0, 1.0};
        double[] K = {3.0, 2.2};

        double[][] Acl = new double[2][2];
        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 2; ++j) {
                Acl[i][j] = A[i][j] - B[i] * K[j];
            }
        }

        // Bryson-style tolerances: x1_max=1, x2_max=2, u_max=0.5.
        double[][] Q = { {1.0, 0.0}, {0.0, 1.0 / 4.0} };
        double R = 1.0 / (0.5 * 0.5);

        if (!isPsd2x2(Q)) {
            throw new IllegalArgumentException("Q must be positive semidefinite.");
        }
        if (R <= 0.0) {
            throw new IllegalArgumentException("R must be positive definite.");
        }

        double[] x = {1.0, 0.0};
        double h = 0.001;
        double tf = 8.0;
        int steps = (int) Math.round(tf / h);
        double J = 0.0;

        double fOld = quad(x, Q) + R * input(K, x) * input(K, x);
        for (int k = 0; k < steps; ++k) {
            double[] xNext = rk4Step(Acl, x, h);
            double uNext = input(K, xNext);
            double fNew = quad(xNext, Q) + R * uNext * uNext;
            J += 0.5 * h * (fOld + fNew);
            x = xNext;
            fOld = fNew;
        }

        System.out.println("Q = [[1, 0], [0, 0.25]]");
        System.out.println("R = " + R);
        System.out.println("Finite-horizon weighted cost J_T = " + J);
    }
}

The MATLAB script uses ode45 for simulation and, when Simulink is available, programmatically creates a small model that computes the state quadratic term \( x^\top Qx \).

Chapter28_Lesson3.m

% Chapter28_Lesson3.m
% Modern Control — Chapter 28, Lesson 3
% State and control weighting matrices Q and R as performance descriptors.
%
% This script demonstrates:
% 1. Bryson-style weight construction.
% 2. Finite-horizon weighted cost for a closed-loop state trajectory.
% 3. Coordinate-transformation consistency.
% 4. Optional programmatic creation of a small Simulink model.

clear; clc;

A = [0 1; -2 -0.4];
B = [0; 1];
K = [3.0 2.2];
x0 = [1; 0];

% Bryson-style tolerances.
xMax = [1; 2];
uMax = 0.5;
Q = diag(1 ./ (xMax.^2));
R = 1 / (uMax^2);

fprintf('eig(Q) = '); disp(eig((Q + Q')/2)');
fprintf('eig(R) = %.6f\n', R);

Acl = A - B*K;
[t, X] = ode45(@(t,x) Acl*x, [0 8], x0);
U = -(K * X')';

stateTerms = sum((X*Q).*X, 2);
inputTerms = (U.^2) * R;
J = trapz(t, stateTerms + inputTerms);
fprintf('Finite-horizon weighted cost J_T = %.6f\n', J);

% Coordinate transformation x = T z. Then Qz = T' Q T.
T = [2 0.5; 0 1.5];
Z = (T \ X')';
Qz = T' * Q * T;
Jx = trapz(t, sum((X*Q).*X, 2));
Jz = trapz(t, sum((Z*Qz).*Z, 2));
fprintf('State cost in x-coordinates = %.6f\n', Jx);
fprintf('State cost in z-coordinates = %.6f\n', Jz);
fprintf('Coordinate consistency error = %.3e\n', abs(Jx - Jz));

% Optional Simulink model: run only when Simulink is available.
if license('test', 'Simulink')
    modelName = 'Chapter28_Lesson3_Simulink_QR_Cost';
    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end
    new_system(modelName);
    open_system(modelName);

    add_block('simulink/Sources/Constant', [modelName '/x_vector'], ...
        'Value', '[1; 0]', 'Position', [80 80 150 120]);
    add_block('simulink/Math Operations/Gain', [modelName '/Q_gain'], ...
        'Gain', mat2str(Q), 'Multiplication', 'Matrix(K*u)', 'Position', [220 70 300 130]);
    add_block('simulink/Math Operations/Dot Product', [modelName '/xTQx'], ...
        'Position', [380 70 460 130]);
    add_block('simulink/Sinks/Display', [modelName '/Display_xTQx'], ...
        'Position', [540 80 650 120]);

    add_line(modelName, 'x_vector/1', 'Q_gain/1');
    add_line(modelName, 'x_vector/1', 'xTQx/1');
    add_line(modelName, 'Q_gain/1', 'xTQx/2');
    add_line(modelName, 'xTQx/1', 'Display_xTQx/1');

    save_system(modelName);
    fprintf('Created optional Simulink model: %s.slx\n', modelName);
end

The Mathematica notebook expression evaluates the same weighted-cost example and symbolically verifies the coordinate-transformation formula.

Chapter28_Lesson3.nb

Notebook[{Cell["Chapter28_Lesson3.nb", "Title"],
Cell["Modern Control — Chapter 28, Lesson 3: State and control weighting matrices Q and R as performance descriptors.", "Text"],
Cell[BoxData[RowBox[{"ClearAll", "[", "\"Global`*\"", "]"}]], "Input"],
Cell[BoxData[RowBox[{RowBox[{"A", "=", RowBox[{"{ {", RowBox[{RowBox[{"0", ",", "1"}], ",", RowBox[{RowBox[{"-", "2"}], ",", RowBox[{"-", "0.4"}]}]}], "} }"}]}], ";"}]], "Input"],
Cell[BoxData[RowBox[{RowBox[{"B", "=", RowBox[{"{ {", RowBox[{"0"}], "},{", RowBox[{"1"}], "} }"}]}], ";", RowBox[{"K", "=", RowBox[{"{ {", RowBox[{"3.0", ",", "2.2"}], "} }"}]}], ";"}]], "Input"],
Cell[BoxData[RowBox[{RowBox[{"xMax", "=", RowBox[{"{", RowBox[{"1", ",", "2"}], "}"}]}], ";", RowBox[{"uMax", "=", "0.5"}], ";"}]], "Input"],
Cell[BoxData[RowBox[{RowBox[{"Q", "=", RowBox[{"DiagonalMatrix", "[", RowBox[{"1", "/", RowBox[{"xMax", "^", "2"}]}], "]"}]}], ";", RowBox[{"R", "=", RowBox[{"{ {", RowBox[{"1", "/", RowBox[{"uMax", "^", "2"}]}], "} }"}]}], ";"}]], "Input"],
Cell[BoxData[RowBox[{"{", RowBox[{RowBox[{"Eigenvalues", "[", "Q", "]"}], ",", RowBox[{"Eigenvalues", "[", "R", "]"}]}], "}"}]], "Input"],
Cell[BoxData[RowBox[{RowBox[{"Acl", "=", RowBox[{"A", "-", RowBox[{"B", ".", "K"}]}]}], ";"}]], "Input"],
Cell[BoxData[RowBox[{"sol", "=", RowBox[{"NDSolveValue", "[", RowBox[{RowBox[{"{", RowBox[{RowBox[{RowBox[{"x", "'"}], "[", "t", "]"}], "==", RowBox[{"Acl", ".", RowBox[{"x", "[", "t", "]"}]}]}], ",", RowBox[{RowBox[{"x", "[", "0", "]"}], "==", RowBox[{"{", RowBox[{"1", ",", "0"}], "}"}]}]}], "}"}], ",", "x", ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "8"}], "}"}]}], "]"}]}]], "Input"],
Cell[BoxData[RowBox[{"J", "=", RowBox[{"NIntegrate", "[", RowBox[{RowBox[{RowBox[{RowBox[{"sol", "[", "t", "]"}], ".", "Q", ".", RowBox[{"sol", "[", "t", "]"}]}], "+", RowBox[{RowBox[{"(", RowBox[{RowBox[{"-", "K"}], ".", RowBox[{"sol", "[", "t", "]"}]}], ")"}], ".", "R", ".", RowBox[{"(", RowBox[{RowBox[{"-", "K"}], ".", RowBox[{"sol", "[", "t", "]"}]}], ")"}]}]}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "8"}], "}"}]}], "]"}]}]], "Input"],
Cell[BoxData[RowBox[{RowBox[{"Tmat", "=", RowBox[{"{ {", RowBox[{RowBox[{"2", ",", "0.5"}], ",", RowBox[{"0", ",", "1.5"}]}], "} }"}]}], ";", RowBox[{"Qz", "=", RowBox[{RowBox[{"Transpose", "[", "Tmat", "]"}], ".", "Q", ".", "Tmat"}]}], ";"}]], "Input"],
Cell[BoxData[RowBox[{"Simplify", "[", RowBox[{RowBox[{RowBox[{"x", ".", "Q", ".", "x"}], "/.", RowBox[{"x", "->", RowBox[{"Tmat", ".", "z"}]}]}], "==", RowBox[{"z", ".", "Qz", ".", "z"}]}], "]"}]], "Input"]
}]

Problem 1 (Symmetric Part): Let \( Q\in\mathbb{R}^{n\times n} \) be any real square matrix. Prove that the quadratic form \( x^\top Qx \) depends only on the symmetric part of \( Q \).

Solution: Define

\[ Q_s=\frac{Q+Q^\top}{2},\qquad Q_a=\frac{Q-Q^\top}{2}. \]

Then \( Q=Q_s+Q_a \) and \( Q_a^\top=-Q_a \). For any real vector \( x \),

\[ x^\top Q_ax=(x^\top Q_ax)^\top=x^\top Q_a^\top x=-x^\top Q_ax. \]

Therefore \( x^\top Q_ax=0 \), so \( x^\top Qx=x^\top Q_sx \).

Problem 2 (Bryson-Style Weights): A state vector is \( x=[x_1,x_2,x_3]^\top \). The acceptable magnitudes are \( x_1^{\max}=2 \), \( x_2^{\max}=0.5 \), and \( x_3^{\max}=10 \). Construct a diagonal \( Q \) using normalized squared penalties.

Solution:

\[ Q=\operatorname{diag}\left(\frac{1}{2^2}, \frac{1}{0.5^2},\frac{1}{10^2}\right) =\operatorname{diag}(0.25,4,0.01). \]

The second state receives the largest weight because its acceptable magnitude is smallest.

Problem 3 (Coordinate Transformation): Suppose \( x=Tz \) and a state cost is \( x^\top Q_xx \). Derive the weighting matrix \( Q_z \) that gives the same scalar cost in \( z \)-coordinates.

Solution:

\[ x^\top Q_xx=(Tz)^\top Q_x(Tz)=z^\top T^\top Q_xTz. \]

Hence \( Q_z=T^\top Q_xT \). This transformed matrix preserves the physical meaning of the original cost.

Problem 4 (Induced Output Weight): Let the performance output be \( z=C_zx \), and suppose the desired output penalty is \( z^\top W_zz \). Find the equivalent state weighting matrix.

Solution:

\[ z^\top W_zz=(C_zx)^\top W_z(C_zx)=x^\top C_z^\top W_zC_zx. \]

Therefore \( Q=C_z^\top W_zC_z \). If \( W_z\succeq 0 \), then \( Q\succeq 0 \) because \( x^\top Qx=(C_zx)^\top W_z(C_zx)\ge 0 \).

Problem 5 (Positive Definiteness of a Two-Input R): Consider \( R=\begin{bmatrix}4&1\\1&2\end{bmatrix} \). Decide whether it is positive definite.

Solution: For a symmetric two-by-two matrix, positive definiteness follows from positive leading principal minors:

\[ 4>0, \qquad \det(R)=4\cdot 2-1\cdot 1=7>0. \]

Thus \( R\succ 0 \), so every nonzero two-input vector has strictly positive input cost.

The matrices \( Q \) and \( R \) are not arbitrary tuning symbols; they are mathematical encodings of state importance, actuator effort, units, physical tolerances, and coordinate choices. A valid quadratic performance descriptor normally uses \( Q\succeq 0 \) and \( R\succ 0 \). Diagonal weights are useful for normalized first designs, while full matrices capture coupled variables or coupled actuators. Coordinate transformations require congruence transformations such as \( Q_z=T^\top Q_xT \), preserving the physical meaning of the quadratic cost. These ideas prepare students for the qualitative weight-behavior discussion in Lesson 4 and the optimal-feedback preview in Lesson 5.