In Chapter 5 we frequently construct a state vector from derivatives, e.g., \( x_1(t) = y(t) \), \( x_2(t) = \dot{y}(t) \), …, \( x_n(t) = y^{(n-1)}(t) \). Even for well-behaved physical systems, these components can have vastly different units and magnitudes (meters vs meters/second vs meters/second²). Such disparities can:

• obscure physical interpretation (states are not comparable on the same numerical scale),
• amplify finite-precision effects in simulation and computation (ill-conditioned operations),
• create implementation issues (fixed-point, ADC/DAC scaling, overflow margins),
• degrade parameter identifiability in regression-based modeling (preview; not pursued here).

The remedy is to define a scaled (often dimensionless) state vector \( z(t) \) related by an invertible linear map \( x(t) = S z(t) \), where \( S \) is typically diagonal.

flowchart TD
  A["Start from ODE or physical model"] --> B["Choose initial state definition (e.g., derivatives)"]
  B --> C["Inspect units and typical magnitudes of each state"]
  C --> D["Pick scaling S (often diagonal): x = S z"]
  D --> E["Transform matrices: A_s = S^{-1} A S, B_s = S^{-1} B, C_s = C S"]
  E --> F["Simulate / compute using scaled coordinates z"]
  F --> G["Map back to physical states: x = S z"]
        

Consider the continuous-time LTI state-space model already introduced in Chapter 4:

\[ \dot{x}(t) = A x(t) + B u(t), \quad y(t) = C x(t) + D u(t), \quad x(t)\in\mathbb{R}^n,\; u(t)\in\mathbb{R}^m,\; y(t)\in\mathbb{R}^p. \]

Let \( S \in \mathbb{R}^{n\times n} \) be nonsingular and define the scaled state \( z(t) \) by \( x(t) = S z(t) \). Differentiating gives \( \dot{x}(t)=S\dot{z}(t) \). Substituting into the state equation:

\[ S\dot{z}(t) = A S z(t) + B u(t) \;\;\Longrightarrow\;\; \dot{z}(t) = \underbrace{S^{-1} A S}_{A_s} z(t) + \underbrace{S^{-1} B}_{B_s} u(t). \]

Likewise, the output becomes:

\[ y(t) = C S z(t) + D u(t) = \underbrace{(C S)}_{C_s} z(t) + D u(t). \]

Thus the scaled realization is \( (A_s,B_s,C_s,D) \) with \( A_s = S^{-1}AS \), \( B_s=S^{-1}B \), \( C_s=CS \).

Theorem 1 (Trajectory Equivalence Under Scaling). Suppose \( S \) is nonsingular. Then for any input \( u(\cdot) \) and any initial condition \( x(0)=x_0 \), the pair of systems \( (A,B,C,D) \) in \( x \)-coordinates and \( (A_s,B_s,C_s,D) \) in \( z \)-coordinates generate the same output trajectory \( y(t) \), provided the initial conditions are related by \( x_0 = S z_0 \).

Proof. Let \( z_0 = S^{-1}x_0 \). If \( z(t) \) solves \( \dot{z} = A_s z + B_s u \) with \( z(0)=z_0 \), define \( x(t)=S z(t) \). Then \( \dot{x}(t)=S\dot{z}(t)=S(A_s z(t)+B_s u(t))=S(S^{-1}AS z(t)+S^{-1}B u(t)) \)
\( =A(S z(t))+B u(t)=A x(t)+B u(t) \) . Therefore \( x(t) \) satisfies the original state equation with \( x(0)=S z_0=x_0 \). . For the output, \( y(t)=C x(t)+D u(t)=C S z(t)+D u(t)=C_s z(t)+D u(t) \), which matches the scaled output equation. ∎

Two important invariants follow immediately:

• Eigenvalues of \(A\) are preserved. Since \( A_s = S^{-1}AS \) is similar to \(A\), they have identical characteristic polynomials and eigenvalues.
• Input-output behavior is unchanged. The mapping \( x = S z \) is internal; outputs match for matched initial conditions. (A formal transfer-function equivalence is deferred to Chapter 6.)

Corollary 1 (Characteristic Polynomial Invariance).

\[ \det(\lambda I - A_s) = \det(\lambda I - S^{-1} A S) = \det(S^{-1}(\lambda I - A)S) = \\ \det(S^{-1})\det(\lambda I - A)\det(S) = \det(\lambda I - A). \]

In engineering practice, \( S \) is often diagonal: \( S=\mathrm{diag}(s_1,\dots,s_n) \) with \( s_i \neq 0 \). Then each state is rescaled independently: \( x_i = s_i z_i \).

A common objective is to select \( s_i \) so that the scaled states satisfy \( |z_i(t)| \) is typically \( \mathcal{O}(1) \) over the operating regime. If a typical (or allowable) magnitude of \( x_i \) is \( X_i \), choose: \( s_i = X_i \), so \( z_i = x_i / X_i \).

Example (Second-Order Mechanical System). Consider the standard mass-spring-damper model:

\[ m\ddot{y}(t) + c\dot{y}(t) + k y(t) = u(t). \]

Using the phase-variable state choice from earlier lessons: \( x_1=y \), \( x_2=\dot{y} \). Then:

\[ \dot{x}(t) = \begin{bmatrix} 0 & 1 \\ -\frac{k}{m} & -\frac{c}{m} \end{bmatrix} x(t) + \begin{bmatrix} 0\\ \frac{1}{m} \end{bmatrix} u(t). \]

Let \( Y \) be a characteristic displacement scale and \( V \) a characteristic velocity scale. Define the dimensionless state \( z \) using \( x = S z \) with \( S=\mathrm{diag}(Y,V) \), i.e., \( z_1=y/Y \), \( z_2=\dot{y}/V \). Then the scaled matrices are:

\[ A_s = S^{-1} A S = \begin{bmatrix} 0 & \frac{V}{Y} \\ -\frac{k}{m}\frac{Y}{V} & -\frac{c}{m} \end{bmatrix}, \quad B_s = S^{-1}B = \begin{bmatrix} 0 \\ \frac{1}{mV} \end{bmatrix}. \]

A particularly informative choice is to set \( V=\omega_0 Y \) where \( \omega_0 = \sqrt{k/m} \) is the undamped natural frequency (known from classical control/mechanics). Then \( V/Y = \omega_0 \) and \( (k/m)(Y/V)=\omega_0 \), yielding a normalized structure:

\[ A_s = \begin{bmatrix} 0 & \omega_0 \\ -\omega_0 & -\frac{c}{m} \end{bmatrix}. \]

If we additionally scale time by \( t = T\, \tau \) with \( T = 1/\omega_0 \), then \( \frac{d}{dt} = \frac{1}{T}\frac{d}{d\tau} = \omega_0 \frac{d}{d\tau} \) and the dynamics become dimensionless in time:

\[ \frac{d z}{d\tau} = \underbrace{\left(\frac{1}{\omega_0}A_s\right)}_{\bar{A}} z + \underbrace{\left(\frac{1}{\omega_0}B_s\right)}_{\bar{B}} u. \]

Time scaling is optional in this course sequence; the key point is that state scaling alone already improves interpretability and computation.

Consider the scalar n-th order LTI ODE (single-input) used throughout Chapter 5:

\[ y^{(n)}(t) + a_{n-1}y^{(n-1)}(t) + \cdots + a_1 \dot{y}(t) + a_0 y(t) = b\,u(t). \]

The unscaled phase-variable choice is \( x_1=y \), \( x_2=\dot{y} \), …, \( x_n = y^{(n-1)} \), producing the companion-form matrix with last row containing \( -a_0,\dots,-a_{n-1} \). When coefficients span orders of magnitude, the resulting \(A\) can be numerically awkward.

A systematic normalization is to introduce a frequency scale \( \omega > 0 \) and define:

\[ z_i(t) = \frac{x_i(t)}{\omega^{i-1}} \quad\Longleftrightarrow\quad x_i(t) = \omega^{i-1} z_i(t), \qquad i=1,\dots,n. \]

This corresponds to a diagonal scaling matrix \( S=\mathrm{diag}(1,\omega,\omega^2,\dots,\omega^{n-1}) \).

Proposition 1 (Structure of the Scaled Companion Matrix). Let \(A\) be the companion matrix for the unscaled phase-variable form. Then \( A_s = S^{-1}AS \) has the same “shift” structure but the last row becomes normalized by powers of \(\omega\).

Proof sketch (direct computation). In the companion form, for \( i=1,\dots,n-1 \) we have \( \dot{x}_i = x_{i+1} \). Using \( x_i = \omega^{i-1} z_i \) gives:

\[ \omega^{i-1}\dot{z}_i = \omega^{i} z_{i+1} \;\;\Longrightarrow\;\; \dot{z}_i = \omega z_{i+1}, \quad i=1,\dots,n-1. \]

Hence the superdiagonal “1” entries become \( \omega \) in the scaled coordinates. For the last equation, \( \dot{x}_n = -a_0 x_1 - a_1 x_2 - \cdots - a_{n-1} x_n + b u \), substituting \( x_k=\omega^{k-1}z_k \) yields:

\[ \omega^{n-1}\dot{z}_n = -\sum_{k=1}^{n} a_{k-1}\,\omega^{k-1} z_k + b u \;\;\Longrightarrow\;\; \dot{z}_n = -\sum_{k=1}^{n} a_{k-1}\,\omega^{k-n} z_k + \frac{b}{\omega^{n-1}}u. \]

Thus coefficients become \( a_{k-1}\omega^{k-n} \), which can be made more uniform by an appropriate choice of \(\omega\). ∎

Choosing \(\omega\). In this chapter, a practical choice is to set \(\omega\) near a dominant system frequency scale inferred from the ODE coefficients (e.g., for a second-order standard form, \(\omega=\sqrt{a_0}\) after proper normalization). The guiding principle is: choose \(\omega\) so that the scaled coefficients have comparable magnitudes.

Although the mathematical system is unchanged (Theorem 1), computations are performed in finite precision. Scaling can reduce roundoff amplification by reducing extreme element magnitudes and improving matrix conditioning. We introduce only the minimal linear-algebra concept needed here.

For a nonsingular matrix \(M\), define its (induced) condition number in a chosen norm: \( \kappa(M)=\|M\|\cdot\|M^{-1}\| \). Large \(\kappa(M)\) indicates that small relative perturbations in data can cause large relative perturbations in results.

Lemma 1 (Similarity Scaling Does Not Worsen Eigenvalues but Can Change Sensitivity). Since \(A_s=S^{-1}AS\), eigenvalues are identical; however, numerical algorithms for eigenvalues and simulation can be more stable on \(A_s\) when \(S\) reduces magnitude disparities.

One simple bound (useful as intuition) relates perturbations through scaling: if \( x = S z \), then relative errors satisfy approximately \( \frac{\|\delta x\|}{\|x\|} \lesssim \kappa(S)\frac{\|\delta z\|}{\|z\|} \). Thus, if we compute in the coordinate system where signals are closer to unity and \(\kappa(S)\) is moderate, we typically obtain improved numerical behavior.

In state-space computations, scaling is especially helpful in:

• time-domain simulation (numerical integration of \(\dot{x}=Ax+Bu\)),
• forming powers/series of \(A\) (used implicitly by some solvers),
• discretization routines (e.g., matrix exponential computations introduced later in Chapter 7–8),
• embedded implementations with limited numeric range.

Without invoking later modern-control concepts, you can select scalings in a disciplined way using only: (i) units, (ii) expected magnitude ranges, and (iii) coefficient normalization from the originating ODE.

Rule A (Unit normalization). Choose \(s_i\) so that each scaled state \(z_i=x_i/s_i\) is dimensionless (or uses consistent base units). This clarifies interpretation and reduces unit-driven magnitude gaps.

Rule B (Range matching). If you have expected maxima \( |x_i(t)| \le X_i \) in operation, set \(s_i=X_i\) so that \( |z_i(t)| \le 1 \) typically.

Rule C (Derivative-based ODE normalization). For phase-variable states derived from derivatives, choose a frequency scale \(\omega\) and define \( x_i=\omega^{i-1}z_i \). This keeps successive derivatives on comparable numerical scales when the signal is dominated by time constants near \(1/\omega\).

In all cases, verify that \(S\) is well-conditioned (avoid extremely large or small \(s_i\) ratios unless required).

We demonstrate scaling on an ODE whose companion realization has coefficients with disparate magnitudes:

\[ y^{(3)}(t) + 1000\,y^{(2)}(t) + 10\,\dot{y}(t) + y(t) = u(t). \]

Define phase-variable states \(x_1=y\), \(x_2=\dot{y}\), \(x_3=y^{(2)}\). Then:

\[ \dot{x}(t) = \begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ -1 & -10 & -1000 \end{bmatrix}x(t) + \begin{bmatrix} 0\\0\\1 \end{bmatrix}u(t),\quad y(t)= \begin{bmatrix} 1 & 0 & 0 \end{bmatrix}x(t). \]

We scale using \(S=\mathrm{diag}(1,\omega,\omega^2)\) with a modest \(\omega\) (e.g., \(\omega=10\)) to reduce derivative magnitude disparities. Libraries relevant to modern control workflows include numpy, scipy, and the python-control package (control).

import numpy as np

# Unscaled companion-form matrices
A = np.array([[0.0, 1.0, 0.0],
              [0.0, 0.0, 1.0],
              [-1.0, -10.0, -1000.0]])
B = np.array([[0.0],
              [0.0],
              [1.0]])
C = np.array([[1.0, 0.0, 0.0]])
D = np.array([[0.0]])

# Choose a frequency scaling omega and build S = diag(1, omega, omega^2)
omega = 10.0
S = np.diag([1.0, omega, omega**2])
Sinv = np.linalg.inv(S)

# Scaled realization
A_s = Sinv @ A @ S
B_s = Sinv @ B
C_s = C @ S
D_s = D.copy()

print("A_s =\n", A_s)
print("B_s =\n", B_s)
print("C_s =\n", C_s)

# Simulate both systems with the same physical initial condition x0 and input u(t)
from scipy.integrate import solve_ivp

def u_of_t(t):
    # unit step input (piecewise constant)
    return 1.0

def f_unscaled(t, x):
    return (A @ x.reshape(-1,1) + B * u_of_t(t)).flatten()

def f_scaled(t, z):
    return (A_s @ z.reshape(-1,1) + B_s * u_of_t(t)).flatten()

t0, tf = 0.0, 0.02
t_eval = np.linspace(t0, tf, 400)

x0 = np.array([0.0, 0.0, 0.0])              # physical initial state
z0 = (Sinv @ x0.reshape(-1,1)).flatten()    # consistent scaled initial state

sol_x = solve_ivp(f_unscaled, (t0, tf), x0, t_eval=t_eval, rtol=1e-9, atol=1e-12)
sol_z = solve_ivp(f_scaled,   (t0, tf), z0, t_eval=t_eval, rtol=1e-9, atol=1e-12)

# Recover x from z: x = S z
x_from_z = (S @ sol_z.y).T
y_unscaled = (C @ sol_x.y).flatten()
y_scaled   = (C @ x_from_z.T).flatten()

# The two outputs should match (up to solver tolerances)
max_err = np.max(np.abs(y_unscaled - y_scaled))
print("max |y_unscaled - y_scaled| =", max_err)

# Optional: use python-control for LTI object handling
try:
    import control
    sys_unscaled = control.ss(A, B, C, D)
    sys_scaled   = control.ss(A_s, B_s, C_s, D_s)
    # Note: internal states differ; outputs coincide when mapped consistently.
    print("Created control.ss objects successfully.")
except Exception as e:
    print("python-control not available in this environment:", e)

In C++, a standard modern linear-algebra choice is the header-only Eigen library. The core operations are: \(A_s=S^{-1}AS\), \(B_s=S^{-1}B\), \(C_s=CS\).

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

int main() {
  using Eigen::MatrixXd;
  using Eigen::VectorXd;

  MatrixXd A(3,3), B(3,1), C(1,3), D(1,1);
  A << 0.0, 1.0, 0.0,
       0.0, 0.0, 1.0,
      -1.0,-10.0,-1000.0;

  B << 0.0, 0.0, 1.0;
  C << 1.0, 0.0, 0.0;
  D << 0.0;

  double omega = 10.0;
  MatrixXd S = MatrixXd::Zero(3,3);
  S(0,0) = 1.0;
  S(1,1) = omega;
  S(2,2) = omega*omega;

  MatrixXd Sinv = S.inverse();

  MatrixXd As = Sinv * A * S;
  MatrixXd Bs = Sinv * B;
  MatrixXd Cs = C * S;
  MatrixXd Ds = D;

  std::cout << "As =\n" << As << "\n\n";
  std::cout << "Bs =\n" << Bs << "\n\n";
  std::cout << "Cs =\n" << Cs << "\n\n";

  // Example: map a scaled state z back to physical x = S z
  VectorXd z(3);
  z << 0.1, -0.2, 0.05;
  VectorXd x = S * z;
  std::cout << "x = S z =\n" << x << "\n";

  return 0;
}

For simulation in C++ you would combine these operations with a numerical ODE integrator (e.g., Runge–Kutta 4 implemented from scratch or libraries such as Boost.Odeint). The mathematical requirement is simply to integrate \(\dot{z}=A_s z + B_s u\) and map back \(x=Sz\).

In Java, EJML is a common numerical linear algebra library. The following code performs the similarity scaling.

import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.CommonOps_DDRM;

public class StateScalingDemo {
  public static void main(String[] args) {
    // A (3x3)
    DMatrixRMaj A = new DMatrixRMaj(new double[][]{
      {0.0, 1.0, 0.0},
      {0.0, 0.0, 1.0},
      {-1.0, -10.0, -1000.0}
    });

    // B (3x1), C (1x3), D (1x1)
    DMatrixRMaj B = new DMatrixRMaj(new double[][]{ {0.0},{0.0},{1.0} });
    DMatrixRMaj C = new DMatrixRMaj(new double[][]{ {1.0,0.0,0.0} });
    DMatrixRMaj D = new DMatrixRMaj(new double[][]{ {0.0} });

    double omega = 10.0;
    DMatrixRMaj S = new DMatrixRMaj(3,3);
    CommonOps_DDRM.setIdentity(S);
    S.set(1,1, omega);
    S.set(2,2, omega*omega);

    DMatrixRMaj Sinv = new DMatrixRMaj(3,3);
    CommonOps_DDRM.invert(S, Sinv);

    // As = Sinv * A * S
    DMatrixRMaj tmp = new DMatrixRMaj(3,3);
    DMatrixRMaj As  = new DMatrixRMaj(3,3);
    CommonOps_DDRM.mult(Sinv, A, tmp);
    CommonOps_DDRM.mult(tmp, S, As);

    // Bs = Sinv * B
    DMatrixRMaj Bs = new DMatrixRMaj(3,1);
    CommonOps_DDRM.mult(Sinv, B, Bs);

    // Cs = C * S
    DMatrixRMaj Cs = new DMatrixRMaj(1,3);
    CommonOps_DDRM.mult(C, S, Cs);

    System.out.println("As =\n" + As);
    System.out.println("Bs =\n" + Bs);
    System.out.println("Cs =\n" + Cs);
    System.out.println("D  =\n" + D);
  }
}

As in C++, simulation requires an ODE solver. A straightforward educational option is to implement RK4 and integrate the scaled system \(\dot{z}=A_s z + B_s u\).

MATLAB’s Control System Toolbox represents state-space systems via ss(A,B,C,D). Scaling can be applied directly using matrix formulas, or via coordinate transformations.

% Unscaled system
A = [0 1 0;
     0 0 1;
    -1 -10 -1000];
B = [0;0;1];
C = [1 0 0];
D = 0;

sys = ss(A,B,C,D);

% Scaling matrix S = diag(1, omega, omega^2)
omega = 10;
S = diag([1, omega, omega^2]);

As = inv(S)*A*S;
Bs = inv(S)*B;
Cs = C*S;
Ds = D;

sys_s = ss(As,Bs,Cs,Ds);

% Verify output equivalence via simulation
t = linspace(0,0.02,400);
u = ones(size(t));     % step input
x0 = [0;0;0];
z0 = inv(S)*x0;

[y1, t1, x1] = lsim(sys,   u, t, x0);
[y2, t2, z2] = lsim(sys_s, u, t, z0);

% Map z back to x and compute output from original C for comparison
x_from_z = (S*z2')';
y2_phys = (C*x_from_z')';

max_err = max(abs(y1 - y2_phys));
disp(["max |y1 - y2_phys| = ", num2str(max_err)]);

Simulink implementation guidance (no images).

• Use a State-Space block for the scaled realization \((A_s,B_s,C_s,D)\).
• Insert a Gain block before the State-Space block to compute \(z(0)=S^{-1}x(0)\) if you initialize from physical states.
• After the State-Space block, if you need physical states, add a Gain block with \(S\) to compute \(x=Sz\).
• Compute outputs either directly from the State-Space block output \(y=C_s z + D u\) or by reconstructing \(x\) and using \(y=Cx+Du\).

Mathematica supports state-space modeling and simulation. We compute scaled matrices and verify output equivalence.

(* Unscaled matrices *)
A = { {0, 1, 0},
     {0, 0, 1},
     {-1, -10, -1000} };
B = { {0}, {0}, {1} };
C = { {1, 0, 0} };
D = { {0} };

omega = 10;
S = DiagonalMatrix[{1, omega, omega^2}];
Sinv = Inverse[S];

As = Sinv.A.S;
Bs = Sinv.B;
Cs = C.S;
Ds = D;

As // MatrixForm
Bs // MatrixForm
Cs // MatrixForm

(* Define input u(t) = 1 and simulate both realizations *)
u[t_] := 1;

t0 = 0; tf = 0.02;
x0 = {0, 0, 0};
z0 = Sinv.x0;

(* ODEs *)
xSol = NDSolveValue[
  {x'[t] == A.x[t] + Flatten[B]*u[t], x[t0] == x0},
  x, {t, t0, tf}
];

zSol = NDSolveValue[
  {z'[t] == As.z[t] + Flatten[Bs]*u[t], z[t0] == z0},
  z, {t, t0, tf}
];

y1[t_] := C.xSol[t] + D*u[t];
y2[t_] := C.(S.zSol[t]) + D*u[t];

(* Compare outputs on a grid *)
grid = Table[t, {t, t0, tf, (tf - t0)/400}];
maxErr = Max[Abs[Table[y1[t] - y2[t], {t, grid}]]];
maxErr

Problem 1 (Derive the Scaled Matrices): Consider the state-space system \( \dot{x}=Ax+Bu \), \( y=Cx+Du \). Let \( x=Sz \) where \(S\) is nonsingular. Derive formulas for \(A_s,B_s,C_s,D_s\) in the \(z\)-coordinates.

Solution: Substitute \(x=Sz\) and \(\dot{x}=S\dot{z}\):

\[ S\dot{z}=ASz+Bu \;\;\Longrightarrow\;\; \dot{z}=S^{-1}AS\,z + S^{-1}B\,u, \quad y=CSz+Du. \]

Hence \(A_s=S^{-1}AS\), \(B_s=S^{-1}B\), \(C_s=CS\), and \(D_s=D\).

Problem 2 (Trajectory Equivalence Proof): Prove that if \(z(t)\) satisfies \( \dot{z}=A_s z + B_s u \) with \(z(0)=S^{-1}x_0\), then \(x(t)=S z(t)\) satisfies \( \dot{x}=Ax+Bu \) with \(x(0)=x_0\).

Solution: Differentiate \(x=Sz\) to obtain \(\dot{x}=S\dot{z}\). Substitute the scaled dynamics:

\[ \dot{x}=S(A_s z + B_s u)=S(S^{-1}AS z + S^{-1}B u)=A(Sz)+Bu=Ax+Bu, \]

and \(x(0)=S z(0)=S(S^{-1}x_0)=x_0\). This completes the proof.

Problem 3 (Scaled Companion Form): For the third-order ODE \( y^{(3)} + a_2 y^{(2)} + a_1 \dot{y} + a_0 y = b u \), form the phase-variable realization with \(x_1=y\), \(x_2=\dot{y}\), \(x_3=y^{(2)}\). Then apply the scaling \(S=\mathrm{diag}(1,\omega,\omega^2)\) and compute \(A_s\), \(B_s\).

Solution: Unscaled:

\[ A= \begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ -a_0 & -a_1 & -a_2 \end{bmatrix}, \quad B= \begin{bmatrix} 0\\0\\b \end{bmatrix}. \]

With \(S=\mathrm{diag}(1,\omega,\omega^2)\), \(S^{-1}=\mathrm{diag}(1,\omega^{-1},\omega^{-2})\). Compute:

\[ A_s=S^{-1}AS= \begin{bmatrix} 0 & \omega & 0\\ 0 & 0 & \omega\\ -a_0\omega^{-2} & -a_1\omega^{-1} & -a_2 \end{bmatrix}, \quad B_s=S^{-1}B= \begin{bmatrix} 0\\0\\b\omega^{-2} \end{bmatrix}. \]

Problem 4 (Design a Unit-Based Scaling): Suppose \(x_1\) is measured in meters and typically satisfies \(|x_1|\le 0.02\), while \(x_2\) is measured in meters/second and typically satisfies \(|x_2|\le 2\). Propose a diagonal scaling \(S\) so that the scaled states satisfy roughly \(|z_i|\le 1\).

Solution: Choose \(S=\mathrm{diag}(0.02,\,2)\) so that \(z_1=x_1/0.02\) and \(z_2=x_2/2\), giving typical magnitudes near unity.

Problem 5 (Similarity and Eigenvalues): Show that \(A\) and \(A_s=S^{-1}AS\) have the same eigenvalues.

Solution: Use determinant invariance under similarity:

\[ \det(\lambda I - A_s) = \det(\lambda I - S^{-1}AS) = \det(S^{-1}(\lambda I - A)S) = \det(\lambda I - A). \]

Since the characteristic polynomials are identical, the eigenvalues (roots) coincide.

Scaling and normalization of state variables are principled coordinate transformations of the form \(x=Sz\) that preserve system behavior while improving unit consistency and numerical robustness. We derived the transformed realization \(A_s=S^{-1}AS\), \(B_s=S^{-1}B\), \(C_s=CS\), proved trajectory/output equivalence, and showed how derivative-based phase-variable states can be normalized using frequency-based diagonal scaling. The implementations demonstrated how to apply scaling in major engineering toolchains.