Pedagogical & Educational Derivations

Routh-Hurwitz Stability Engine

routh

Routh-Hurwitz stability criterion and step-by-step array construction.

RouthResult module-attribute

RouthResult = RouthArray

RouthArray

Step-by-step pedagogical Routh-Hurwitz stability analyzer.

Constructs the complete Routh array analytically using SymPy, handles edge cases (leading zero substitution with \(\epsilon > 0\), rows of all zeros with auxiliary polynomial differentiation), evaluates RHP pole counts, and solves for parametric stability intervals.

Parameters:

Name	Type	Description	Default
poly_or_sys	TransferFunction | Sequence[float | Basic] | NDArray | Expr | Poly	Characteristic polynomial coefficients (highest degree first), polynomial expression, or a continuous-time TransferFunction.	required
variable	str | Symbol	Laplace complex frequency variable symbol.	"s"
k_symbol	str | Symbol | None	Optional parametric gain symbol to solve for closed-loop stability intervals.	None

Attributes:

Name	Type	Description
table	list[list[Any]]	The complete computed 2D Routh array containing numeric or symbolic expressions.
row_labels	list[str]	Labels for each row of the Routh array (e.g. ['s^3', 's^2', 's^1', 's^0']).
num_rhp_poles	int	Number of roots in the open Right-Half Plane (RHP). Equal to sign changes in column 1.
is_stable	bool	True if all roots lie strictly in the open Left-Half Plane (LHP).
k_range	str | None	Parametric stability range for symbolic gain \(K\), if applicable.
auxiliary_polynomials	list[Expr]	Any auxiliary polynomials formed when encountering rows of all zeros.
steps	list[str]	Pedagogical step-by-step construction notes.
polynomial	Any	The original characteristic polynomial analyzed.
variable	Symbol	The symbolic variable used for polynomial powers.

Source code in src/ctrlpy/symbolic/routh.py

class RouthArray:
    r"""Step-by-step pedagogical Routh-Hurwitz stability analyzer.

    Constructs the complete Routh array analytically using SymPy, handles edge cases
    (leading zero substitution with $\epsilon > 0$, rows of all zeros with auxiliary polynomial
    differentiation), evaluates RHP pole counts, and solves for parametric stability intervals.

    Parameters
    ----------
    poly_or_sys : TransferFunction | Sequence[float | sp.Basic] | NDArray | sp.Expr | sp.Poly
        Characteristic polynomial coefficients (highest degree first), polynomial expression,
        or a continuous-time TransferFunction.
    variable : str | sp.Symbol, default="s"
        Laplace complex frequency variable symbol.
    k_symbol : str | sp.Symbol | None, default=None
        Optional parametric gain symbol to solve for closed-loop stability intervals.

    Attributes
    ----------
    table : list[list[Any]]
        The complete computed 2D Routh array containing numeric or symbolic expressions.
    row_labels : list[str]
        Labels for each row of the Routh array (e.g. `['s^3', 's^2', 's^1', 's^0']`).
    num_rhp_poles : int
        Number of roots in the open Right-Half Plane (RHP). Equal to sign changes in column 1.
    is_stable : bool
        True if all roots lie strictly in the open Left-Half Plane (LHP).
    k_range : str | None
        Parametric stability range for symbolic gain $K$, if applicable.
    auxiliary_polynomials : list[sp.Expr]
        Any auxiliary polynomials formed when encountering rows of all zeros.
    steps : list[str]
        Pedagogical step-by-step construction notes.
    polynomial : Any
        The original characteristic polynomial analyzed.
    variable : sp.Symbol
        The symbolic variable used for polynomial powers.
    """

    def __init__(
        self,
        poly_or_sys: Any,
        variable: str | sp.Symbol = "s",
        k_symbol: str | sp.Symbol | None = None,
    ) -> None:
        s = sp.Symbol(variable) if isinstance(variable, str) else variable
        k_sym = sp.Symbol(k_symbol) if isinstance(k_symbol, str) else k_symbol

        steps: list[str] = []
        detailed_derivations: list[str] = []
        auxiliary_polys: list[Any] = []

        # 1. Parse input into list of coefficients
        if isinstance(poly_or_sys, TransferFunction):
            if k_sym is not None:
                # Closed loop characteristic polynomial: D(s) + K * N(s)
                num_poly = sum(
                    c * s ** (len(poly_or_sys.num) - 1 - idx)
                    for idx, c in enumerate(poly_or_sys.num)
                )
                den_poly = sum(
                    c * s ** (len(poly_or_sys.den) - 1 - idx)
                    for idx, c in enumerate(poly_or_sys.den)
                )
                char_expr = sp.expand(den_poly + k_sym * num_poly)
                p = sp.Poly(char_expr, s)
                coeffs = [sp.simplify(c) for c in p.all_coeffs()]
                poly_display = char_expr
            else:
                coeffs = [sp.sympify(c) for c in poly_or_sys.den]
                poly_display = sum(c * s ** (len(coeffs) - 1 - idx) for idx, c in enumerate(coeffs))
        elif isinstance(poly_or_sys, (list, tuple, np.ndarray)):
            coeffs = [sp.sympify(c) for c in poly_or_sys]
            poly_display = sum(c * s ** (len(coeffs) - 1 - idx) for idx, c in enumerate(coeffs))
        elif isinstance(poly_or_sys, sp.Poly):
            coeffs = [sp.simplify(c) for c in poly_or_sys.all_coeffs()]
            poly_display = poly_or_sys.as_expr()
        elif isinstance(poly_or_sys, sp.Basic):
            p = sp.Poly(poly_or_sys, s)
            coeffs = [sp.simplify(c) for c in p.all_coeffs()]
            poly_display = poly_or_sys
        else:
            raise TypeError(f"Unsupported polynomial or system type: {type(poly_or_sys)}")

        # Auto-detect parametric gain symbol if present
        if k_sym is None:
            free_syms = set()
            for c in coeffs:
                if isinstance(c, sp.Basic):
                    free_syms.update(c.free_symbols)
            free_syms.discard(s)
            if len(free_syms) == 1:
                k_sym = next(iter(free_syms))

        n = len(coeffs) - 1
        if n < 0:
            raise ValueError("Empty polynomial provided.")

        self.polynomial = poly_display
        self.variable = s
        self.k_symbol = k_sym

        detailed_derivations.append(
            f"**Characteristic Equation:** $P({s}) = {sp.latex(poly_display)} = 0$ (Degree $n = {n}$)"
        )

        if n == 0:
            self.table = [[coeffs[0]]]
            self.num_rhp_poles = 0
            self.is_stable = True
            self.k_range = None
            self.steps = ["Degree 0 polynomial (constant)."]
            self.row_labels = ["s^0"]
            self.auxiliary_polynomials = []
            self._detailed_derivations = detailed_derivations
            return

        num_cols = (n + 2) // 2
        row_labels = [f"s^{n - i}" for i in range(n + 1)]
        table: list[list[Any]] = []

        # Row 0 (s^n) and Row 1 (s^(n-1))
        row0: list[Any] = [
            coeffs[2 * j] if 2 * j < len(coeffs) else sp.Integer(0) for j in range(num_cols)
        ]
        row1: list[Any] = [
            coeffs[2 * j + 1] if 2 * j + 1 < len(coeffs) else sp.Integer(0) for j in range(num_cols)
        ]
        table.append(row0)
        table.append(row1)

        steps.append(
            f"Initialized Row {row_labels[0]} and Row {row_labels[1]} from polynomial coefficients."
        )
        row0_str = ", ".join(f"${sp.latex(x)}$" for x in row0)
        row1_str = ", ".join(f"${sp.latex(x)}$" for x in row1)
        detailed_derivations.append(
            f"**Row ${row_labels[0]}$ Initialization:** [{row0_str}] (alternating even/odd coefficients)"
        )
        detailed_derivations.append(
            f"**Row ${row_labels[1]}$ Initialization:** [{row1_str}] (alternating coefficients)"
        )

        eps = sp.Symbol("epsilon", positive=True)

        # Compute rows 2 through n
        for i in range(2, n + 1):
            prev_row = table[i - 1]
            prev_prev_row = table[i - 2]
            pwr_label = row_labels[i]

            # Special Case 1: First element in previous row is 0, but row is not all zeros
            if _is_zero(prev_row[0]) and not all(_is_zero(x) for x in prev_row):
                table[i - 1][0] = eps
                prev_row[0] = eps
                note1 = (
                    f"Special Case 1: First element in row {row_labels[i - 1]} is zero. "
                    r"Replaced with symbolic $\epsilon > 0$."
                )
                steps.append(note1)
                detailed_derivations.append(f"**{note1}**")

            # Special Case 2: Previous row is all zeros
            if all(_is_zero(x) for x in prev_row):
                aux_pwr = n - (i - 2)
                aux_terms = []
                deriv_row: list[Any] = []
                for j in range(num_cols):
                    term_pwr = aux_pwr - 2 * j
                    c = prev_prev_row[j]
                    if term_pwr >= 0 and not _is_zero(c):
                        aux_terms.append(c * s**term_pwr)
                    deriv_c = c * term_pwr if term_pwr >= 0 else sp.Integer(0)
                    deriv_row.append(sp.simplify(deriv_c))

                aux_poly = sum(aux_terms)
                auxiliary_polys.append(aux_poly)
                table[i - 1] = deriv_row
                prev_row = deriv_row
                note2 = (
                    f"Special Case 2: Row {row_labels[i - 1]} was all zeros. "
                    f"Formed auxiliary polynomial $A({s}) = {sp.latex(aux_poly)}$ from row {row_labels[i - 2]} "
                    r"and replaced zero row with $\frac{dA}{ds}$ coefficients."
                )
                steps.append(note2)
                detailed_derivations.append(f"**{note2}**")

            # Compute row i entries: R[i, j] = (R[i-1, 0]*R[i-2, j+1] - R[i-2, 0]*R[i-1, j+1]) / R[i-1, 0]
            new_row: list[Any] = []
            a = prev_row[0]
            c = prev_prev_row[0]
            row_derivations: list[str] = []

            for j in range(num_cols):
                b = prev_prev_row[j + 1] if j + 1 < num_cols else sp.Integer(0)
                d = prev_row[j + 1] if j + 1 < num_cols else sp.Integer(0)
                if _is_zero(a):
                    entry = sp.Integer(0)
                else:
                    entry = sp.simplify((a * b - c * d) / a)
                new_row.append(entry)
                row_derivations.append(
                    rf"c_{{{j}}} = \frac{{({sp.latex(a)})({sp.latex(b)}) - ({sp.latex(c)})({sp.latex(d)})}}{{{sp.latex(a)}}} = {sp.latex(entry)}"
                )

            # If the newly computed row is all zeros and not at the end
            if all(_is_zero(x) for x in new_row) and i <= n:
                aux_pwr = n - (i - 1)
                aux_terms = []
                deriv_row = []
                for j in range(num_cols):
                    term_pwr = aux_pwr - 2 * j
                    coeff_val = prev_row[j]
                    if term_pwr >= 0 and not _is_zero(coeff_val):
                        aux_terms.append(coeff_val * s**term_pwr)
                    deriv_c = coeff_val * term_pwr if term_pwr >= 0 else sp.Integer(0)
                    deriv_row.append(sp.simplify(deriv_c))

                aux_poly = sum(aux_terms)
                auxiliary_polys.append(aux_poly)
                note2_new = (
                    f"Special Case 2: Row {pwr_label} computed as all zeros. "
                    f"Formed auxiliary polynomial $A({s}) = {sp.latex(aux_poly)}$ from row {row_labels[i - 1]} "
                    r"and replaced with $\frac{dA}{ds}$ coefficients."
                )
                steps.append(note2_new)
                detailed_derivations.append(f"**{note2_new}**")
                new_row = deriv_row

            table.append(new_row)
            entries_str = ", ".join(f"${sp.latex(x)}$" for x in new_row)
            detailed_derivations.append(f"**Row ${pwr_label}$ Computation:** [{entries_str}]")

        # 3. Analyze Column 1 for Stability & Sign Changes
        col1 = [row[0] for row in table]
        k_range = _solve_k_range(col1, k_sym) if k_sym is not None else None

        signs: list[int] = []
        has_unresolved_param = False
        for entry in col1:
            if isinstance(entry, sp.Basic) and k_sym is not None and k_sym in entry.free_symbols:
                has_unresolved_param = True
                break
            sign = _eval_sign_at_limit(entry, eps)
            signs.append(sign)

        if has_unresolved_param:
            num_rhp = -1
            is_stable = False
            steps.append(f"Parametric analysis performed for symbolic gain {k_sym}.")
            detailed_derivations.append(
                rf"**Parametric Gain Analysis:** Solved column 1 positivity for ${k_sym} \implies$ **{k_range}**"
            )
        else:
            num_rhp = 0
            non_zero_signs = [sg for sg in signs if sg != 0]
            for idx in range(len(non_zero_signs) - 1):
                if non_zero_signs[idx] * non_zero_signs[idx + 1] < 0:
                    num_rhp += 1

            is_stable = (
                (num_rhp == 0)
                and (len(auxiliary_polys) == 0)
                and (len(non_zero_signs) == len(signs))
            )
            col1_tex = ", ".join(f"${sp.latex(x)}$" for x in col1)
            steps.append(f"Column 1 signs: {signs} -> {num_rhp} sign changes detected.")
            detailed_derivations.append(
                rf"**First Column Entries:** [{col1_tex}] with signs {signs} $\implies$ **{num_rhp} sign change(s)**."
            )
            conclusion_text = (
                "**Stability Conclusion:** Strictly Asymptotically Stable (All roots in open LHP)."
                if is_stable
                else (
                    f"**Stability Conclusion:** Unstable / Marginally Stable with **{num_rhp} Right-Half Plane (RHP) pole(s)**."
                )
            )
            detailed_derivations.append(conclusion_text)

        self.table = table
        self.num_rhp_poles = num_rhp
        self.is_stable = is_stable
        self.k_range = k_range
        self.steps = steps
        self.row_labels = row_labels
        self.auxiliary_polynomials = auxiliary_polys
        self._detailed_derivations = detailed_derivations

    def explain_steps(self) -> list[str]:
        """Return a step-by-step list of Markdown/LaTeX explanations detailing the array derivation.

        Returns
        -------
        list[str]
            List of explanatory step descriptions for each phase of the
            Routh-Hurwitz construction and stability evaluation.
        """
        return list(self._detailed_derivations)

    def _repr_latex_(self) -> str:
        """Render the Routh table and stability results as formatted LaTeX for Jupyter."""
        max_cols = max((len(r) for r in self.table), default=1)
        col_spec = "c|" + "c" * max_cols
        rows_latex: list[str] = []

        for label, row in zip(self.row_labels, self.table, strict=False):
            lbl_tex = f"{label}"
            if "^" in label:
                base, pwr = label.split("^", 1)
                lbl_tex = f"{base}^{{{pwr}}}"

            entries = [
                sp.latex(val) if hasattr(val, "as_expr") or isinstance(val, sp.Basic) else str(val)
                for val in row
            ]
            rows_latex.append(f"{lbl_tex} & " + " & ".join(entries) + r" \\")

        table_body = "\n".join(rows_latex)
        latex_str = (
            r"\begin{aligned}"
            r"\textbf{Routh-Hurwitz Array:}" + "\n"
            r"\begin{array}{" + col_spec + r"}" + "\n" + table_body + "\n"
            r"\end{array}" + "\n\n"
        )

        status_str = (
            r"\text{Strictly Stable (All LHP poles)}"
            if self.is_stable
            else r"\text{Unstable / Marginally Stable}"
        )
        latex_str += rf"\textbf{{Stability Status:}} &\quad {status_str} \\" + "\n"
        latex_str += rf"\textbf{{RHP Poles:}} &\quad {self.num_rhp_poles} \\" + "\n"

        if self.k_range:
            latex_str += rf"\textbf{{Stable Gain Range (K):}} &\quad {self.k_range} \\" + "\n"

        if self.auxiliary_polynomials:
            aux_str = ", ".join(sp.latex(p) for p in self.auxiliary_polynomials)
            latex_str += rf"\textbf{{Auxiliary Polynomials:}} &\quad {aux_str} \\" + "\n"

        latex_str += r"\end{aligned}"
        return latex_str

    def _repr_markdown_(self) -> str:
        """Return a Markdown representation for Jupyter environments."""
        return f"$${self._repr_latex_()}$$"

    def __str__(self) -> str:
        """Format a human-readable ASCII representation of the Routh table."""
        lines = ["=== Routh-Hurwitz Stability Criterion ==="]
        if self.polynomial is not None:
            lines.append(f"Polynomial: {self.polynomial} = 0")
        lines.append("")

        col_strs: list[list[str]] = []
        for label, row in zip(self.row_labels, self.table, strict=False):
            row_str = [label] + [str(sp.simplify(v) if isinstance(v, sp.Basic) else v) for v in row]
            col_strs.append(row_str)

        max_cols = max(len(r) for r in col_strs)
        col_widths = [0] * max_cols
        for r in col_strs:
            for c_idx, val in enumerate(r):
                col_widths[c_idx] = max(col_widths[c_idx], len(val))

        for r in col_strs:
            label_part = r[0].rjust(col_widths[0])
            entries_part = " | " + "  ".join(r[i].rjust(col_widths[i]) for i in range(1, len(r)))
            lines.append(f"{label_part}{entries_part}")

        lines.append("-" * 42)
        lines.append(f"Number of RHP Poles: {self.num_rhp_poles}")
        lines.append(f"Asymptotically Stable: {self.is_stable}")
        if self.k_range:
            lines.append(f"Stability Range for K: {self.k_range}")
        if self.auxiliary_polynomials:
            lines.append(f"Auxiliary Polynomials: {self.auxiliary_polynomials}")
        if self.steps:
            lines.append("\nConstruction Notes:")
            for s_note in self.steps:
                lines.append(f"  - {s_note}")
        return "\n".join(lines)

explain_steps

explain_steps()

Return a step-by-step list of Markdown/LaTeX explanations detailing the array derivation.

Returns:

Type	Description
list[str]	List of explanatory step descriptions for each phase of the Routh-Hurwitz construction and stability evaluation.

Source code in src/ctrlpy/symbolic/routh.py

def explain_steps(self) -> list[str]:
    """Return a step-by-step list of Markdown/LaTeX explanations detailing the array derivation.

    Returns
    -------
    list[str]
        List of explanatory step descriptions for each phase of the
        Routh-Hurwitz construction and stability evaluation.
    """
    return list(self._detailed_derivations)

__str__

__str__()

Format a human-readable ASCII representation of the Routh table.

Source code in src/ctrlpy/symbolic/routh.py

def __str__(self) -> str:
    """Format a human-readable ASCII representation of the Routh table."""
    lines = ["=== Routh-Hurwitz Stability Criterion ==="]
    if self.polynomial is not None:
        lines.append(f"Polynomial: {self.polynomial} = 0")
    lines.append("")

    col_strs: list[list[str]] = []
    for label, row in zip(self.row_labels, self.table, strict=False):
        row_str = [label] + [str(sp.simplify(v) if isinstance(v, sp.Basic) else v) for v in row]
        col_strs.append(row_str)

    max_cols = max(len(r) for r in col_strs)
    col_widths = [0] * max_cols
    for r in col_strs:
        for c_idx, val in enumerate(r):
            col_widths[c_idx] = max(col_widths[c_idx], len(val))

    for r in col_strs:
        label_part = r[0].rjust(col_widths[0])
        entries_part = " | " + "  ".join(r[i].rjust(col_widths[i]) for i in range(1, len(r)))
        lines.append(f"{label_part}{entries_part}")

    lines.append("-" * 42)
    lines.append(f"Number of RHP Poles: {self.num_rhp_poles}")
    lines.append(f"Asymptotically Stable: {self.is_stable}")
    if self.k_range:
        lines.append(f"Stability Range for K: {self.k_range}")
    if self.auxiliary_polynomials:
        lines.append(f"Auxiliary Polynomials: {self.auxiliary_polynomials}")
    if self.steps:
        lines.append("\nConstruction Notes:")
        for s_note in self.steps:
            lines.append(f"  - {s_note}")
    return "\n".join(lines)

routh_array

routh_array(poly_or_sys, variable='s', k_symbol=None)

Construct a RouthArray analyzer for a polynomial or dynamic system.

Parameters:

Name	Type	Description	Default
poly_or_sys	TransferFunction | Sequence[float | Basic] | NDArray | Expr | Poly	Characteristic polynomial coefficients, expression, or TransferFunction.	required
variable	str | Symbol	Laplace frequency variable symbol.	"s"
k_symbol	str | Symbol | None	Optional parametric gain symbol to solve for stability boundaries.	None

Returns:

Type	Description
RouthArray	The computed Routh-Hurwitz stability array.

Source code in src/ctrlpy/symbolic/routh.py

def routh_array(
    poly_or_sys: Any,
    variable: str | sp.Symbol = "s",
    k_symbol: str | sp.Symbol | None = None,
) -> RouthArray:
    """Construct a RouthArray analyzer for a polynomial or dynamic system.

    Parameters
    ----------
    poly_or_sys : TransferFunction | Sequence[float | sp.Basic] | NDArray | sp.Expr | sp.Poly
        Characteristic polynomial coefficients, expression, or TransferFunction.
    variable : str | sp.Symbol, default="s"
        Laplace frequency variable symbol.
    k_symbol : str | sp.Symbol | None, default=None
        Optional parametric gain symbol to solve for stability boundaries.

    Returns
    -------
    RouthArray
        The computed Routh-Hurwitz stability array.
    """
    return RouthArray(poly_or_sys=poly_or_sys, variable=variable, k_symbol=k_symbol)

routh_table

routh_table(poly_or_sys, variable='s', k_symbol=None)

Construct the Routh-Hurwitz array and analyze system stability.

Parameters:

Name	Type	Description	Default
poly_or_sys	TransferFunction | Sequence[float | Basic] | NDArray | Expr | Poly	The polynomial coefficients (highest degree first), characteristic expression, or a TransferFunction model.	required
variable	str | Symbol	The Laplace variable symbol.	"s"
k_symbol	str | Symbol | None	Optional parametric gain symbol to solve for stability boundaries.	None

Returns:

Type	Description
RouthArray	Computed Routh array object with table, stability metrics, and derivations.

Source code in src/ctrlpy/symbolic/routh.py

def routh_table(
    poly_or_sys: Any,
    variable: str | sp.Symbol = "s",
    k_symbol: str | sp.Symbol | None = None,
) -> RouthArray:
    """Construct the Routh-Hurwitz array and analyze system stability.

    Parameters
    ----------
    poly_or_sys : TransferFunction | Sequence[float | sp.Basic] | NDArray | sp.Expr | sp.Poly
        The polynomial coefficients (highest degree first), characteristic expression,
        or a TransferFunction model.
    variable : str | sp.Symbol, default="s"
        The Laplace variable symbol.
    k_symbol : str | sp.Symbol | None, default=None
        Optional parametric gain symbol to solve for stability boundaries.

    Returns
    -------
    RouthArray
        Computed Routh array object with table, stability metrics, and derivations.
    """
    return RouthArray(poly_or_sys=poly_or_sys, variable=variable, k_symbol=k_symbol)

Analytical Root Locus Rules

root_locus

Classroom analytical Root Locus rules and step-by-step derivations.

RootLocusRulesResult module-attribute

RootLocusRulesResult = RootLocusRules

RootLocusRules

Classroom pedagogical Root Locus rule derivation and analyzer.

Performs complete analytical calculations for all classical Evans root locus rules: 1. Branches and Terminations: Number of open-loop poles (\(n\)), zeros (\(m\)), branches (\(n\)), and branches terminating at infinity (\(n - m\)). 2. Real-Axis Segments: Identifies segments on the real axis lying to the left of an odd number of real poles and zeros. 3. Asymptotes: Calculates the real-axis centroid \(\sigma_a = \frac{\sum p_i - \sum z_j}{n - m}\) and asymptote angles \(\theta_a = \frac{(2k+1)\cdot 180^\circ}{n - m}\). 4. Breakaway and Break-in Points: Solves \(\frac{dK}{ds} = 0 \implies D'(s)N(s) - D(s)N'(s) = 0\) for real points on the locus with \(K \ge 0\). 5. Departure and Arrival Angles: Calculates departure angles from complex open-loop poles \(\theta_d = 180^\circ - \sum \theta_p + \sum \theta_z\) and arrival angles at complex zeros \(\theta_a = 180^\circ - \sum \theta_z + \sum \theta_p\). 6. Imaginary-Axis Crossings: Solves for crossing frequencies \(\omega\) and critical stability gains \(K_{\text{crit}}\) where branches cross the \(j\omega\)-axis (\(s = j\omega\)).

Parameters:

Name	Type	Description	Default
sys_or_poly	TransferFunction | Any	Open-loop transfer function \(G(s)H(s)\), LTI model, or symbolic expression.	None
num	Sequence[float] | None	Numerator coefficients if specified explicitly.	None
den	Sequence[float] | None	Denominator coefficients if specified explicitly.	None
poles	Sequence[complex | float] | None	Explicit list of open-loop poles.	None
zeros	Sequence[complex | float] | None	Explicit list of open-loop zeros.	None
gain	float	Open-loop scalar gain factor.	1.0

Attributes:

Name	Type	Description
num_poles	int	Number of open-loop poles (\(n\)).
num_zeros	int	Number of open-loop zeros (\(m\)).
num_branches	int	Number of Root Locus branches (\(\max(n, m)\)).
num_asymptotes	int	Number of branches approaching infinity (\(\max(0, n - m)\)).
poles	list[complex]	Open-loop pole locations.
zeros	list[complex]	Open-loop zero locations.
real_axis_segments	list[tuple[float, float]]	Intervals on the real axis belonging to the root locus.
centroid	float | None	Asymptote center of mass (\(\sigma_a\)) on the real axis.
asymptote_angles_deg	list[float]	Asymptote angles in degrees.
breakaway_points	list[dict[str, Any]]	Breakaway and break-in points with their critical gain \(K\) and type.
departure_angles	dict[complex, float]	Angle of departure from complex open-loop poles in degrees.
arrival_angles	dict[complex, float]	Angle of arrival to complex open-loop zeros in degrees.
imag_axis_crossings	list[dict[str, Any]]	Crossing frequencies (\(\omega\)) and critical gains (\(K_{\text{crit}}\)) on the imaginary axis.
steps	list[str]	Pedagogical step-by-step summary for each rule.

Source code in src/ctrlpy/symbolic/root_locus.py

class RootLocusRules:
    r"""Classroom pedagogical Root Locus rule derivation and analyzer.

    Performs complete analytical calculations for all classical Evans root locus rules:
    1. **Branches and Terminations:** Number of open-loop poles ($n$), zeros ($m$), branches ($n$),
       and branches terminating at infinity ($n - m$).
    2. **Real-Axis Segments:** Identifies segments on the real axis lying to the left of an odd
       number of real poles and zeros.
    3. **Asymptotes:** Calculates the real-axis centroid $\sigma_a = \frac{\sum p_i - \sum z_j}{n - m}$
       and asymptote angles $\theta_a = \frac{(2k+1)\cdot 180^\circ}{n - m}$.
    4. **Breakaway and Break-in Points:** Solves $\frac{dK}{ds} = 0 \implies D'(s)N(s) - D(s)N'(s) = 0$
       for real points on the locus with $K \ge 0$.
    5. **Departure and Arrival Angles:** Calculates departure angles from complex open-loop poles
       $\theta_d = 180^\circ - \sum \theta_p + \sum \theta_z$ and arrival angles at complex zeros
       $\theta_a = 180^\circ - \sum \theta_z + \sum \theta_p$.
    6. **Imaginary-Axis Crossings:** Solves for crossing frequencies $\omega$ and critical stability
       gains $K_{\text{crit}}$ where branches cross the $j\omega$-axis ($s = j\omega$).

    Parameters
    ----------
    sys_or_poly : TransferFunction | Any, optional
        Open-loop transfer function $G(s)H(s)$, LTI model, or symbolic expression.
    num : Sequence[float] | None, default=None
        Numerator coefficients if specified explicitly.
    den : Sequence[float] | None, default=None
        Denominator coefficients if specified explicitly.
    poles : Sequence[complex | float] | None, default=None
        Explicit list of open-loop poles.
    zeros : Sequence[complex | float] | None, default=None
        Explicit list of open-loop zeros.
    gain : float, default=1.0
        Open-loop scalar gain factor.

    Attributes
    ----------
    num_poles : int
        Number of open-loop poles ($n$).
    num_zeros : int
        Number of open-loop zeros ($m$).
    num_branches : int
        Number of Root Locus branches ($\max(n, m)$).
    num_asymptotes : int
        Number of branches approaching infinity ($\max(0, n - m)$).
    poles : list[complex]
        Open-loop pole locations.
    zeros : list[complex]
        Open-loop zero locations.
    real_axis_segments : list[tuple[float, float]]
        Intervals on the real axis belonging to the root locus.
    centroid : float | None
        Asymptote center of mass ($\sigma_a$) on the real axis.
    asymptote_angles_deg : list[float]
        Asymptote angles in degrees.
    breakaway_points : list[dict[str, Any]]
        Breakaway and break-in points with their critical gain $K$ and type.
    departure_angles : dict[complex, float]
        Angle of departure from complex open-loop poles in degrees.
    arrival_angles : dict[complex, float]
        Angle of arrival to complex open-loop zeros in degrees.
    imag_axis_crossings : list[dict[str, Any]]
        Crossing frequencies ($\omega$) and critical gains ($K_{\text{crit}}$) on the imaginary axis.
    steps : list[str]
        Pedagogical step-by-step summary for each rule.
    """

    def __init__(
        self,
        sys_or_poly: Any = None,
        *,
        num: Sequence[float] | None = None,
        den: Sequence[float] | None = None,
        poles: Sequence[complex | float] | None = None,
        zeros: Sequence[complex | float] | None = None,
        gain: float = 1.0,
    ) -> None:
        steps: list[str] = []
        detailed_derivations: list[str] = []

        # Parse inputs
        parsed_num: np.ndarray | None = None
        parsed_den: np.ndarray | None = None
        parsed_poles: list[complex] = []
        parsed_zeros: list[complex] = []

        if isinstance(sys_or_poly, TransferFunction):
            parsed_num = np.asarray(sys_or_poly.num, dtype=float)
            parsed_den = np.asarray(sys_or_poly.den, dtype=float)
            parsed_poles = [complex(p) for p in sys_or_poly.poles()]
            parsed_zeros = [complex(z) for z in sys_or_poly.zeros()]
        elif sys_or_poly is not None and hasattr(sys_or_poly, "to_tf"):
            tf_obj = sys_or_poly.to_tf()
            parsed_num = np.asarray(tf_obj.num, dtype=float)
            parsed_den = np.asarray(tf_obj.den, dtype=float)
            parsed_poles = [complex(p) for p in tf_obj.poles()]
            parsed_zeros = [complex(z) for z in tf_obj.zeros()]
        elif sys_or_poly is not None and isinstance(sys_or_poly, sp.Basic):
            # Parse symbolic expression
            s = sp.Symbol("s")
            # Remove any symbolic gain K if present in numerator
            expr = sys_or_poly
            numer, denom = sp.fraction(sp.together(expr))
            p_num = sp.Poly(numer, s)
            p_den = sp.Poly(denom, s)
            c_num = [float(c) for c in p_num.all_coeffs()]
            c_den = [float(c) for c in p_den.all_coeffs()]
            parsed_num = np.asarray(c_num, dtype=float)
            parsed_den = np.asarray(c_den, dtype=float)
            parsed_poles = [complex(r) for r in np.roots(parsed_den)] if len(parsed_den) > 1 else []
            parsed_zeros = [complex(r) for r in np.roots(parsed_num)] if len(parsed_num) > 1 else []
        elif num is not None and den is not None:
            parsed_num = np.asarray(num, dtype=float)
            parsed_den = np.asarray(den, dtype=float)
            parsed_poles = [complex(r) for r in np.roots(parsed_den)] if len(parsed_den) > 1 else []
            parsed_zeros = [complex(r) for r in np.roots(parsed_num)] if len(parsed_num) > 1 else []
        elif poles is not None:
            parsed_poles = [complex(p) for p in poles]
            parsed_zeros = [complex(z) for z in zeros] if zeros is not None else []
            # Form polynomials
            p_poly = np.poly(parsed_poles) if parsed_poles else np.array([1.0])
            z_poly = gain * np.poly(parsed_zeros) if parsed_zeros else np.array([float(gain)])
            parsed_num = np.asarray(z_poly, dtype=float)
            parsed_den = np.asarray(p_poly, dtype=float)
        elif isinstance(sys_or_poly, (list, tuple)) and len(sys_or_poly) == 2:
            parsed_num = np.asarray(sys_or_poly[0], dtype=float)
            parsed_den = np.asarray(sys_or_poly[1], dtype=float)
            parsed_poles = [complex(r) for r in np.roots(parsed_den)] if len(parsed_den) > 1 else []
            parsed_zeros = [complex(r) for r in np.roots(parsed_num)] if len(parsed_num) > 1 else []
        else:
            raise TypeError(
                "RootLocusRules expects a TransferFunction, LTI object, (num, den) pair, "
                "or explicit poles/zeros."
            )

        if parsed_num is None or parsed_den is None:
            raise ValueError("Could not construct open-loop numerator and denominator polynomials.")

        # Ensure monic or float representation
        n = len(parsed_poles)
        m = len(parsed_zeros)
        num_branches = max(n, m)
        num_asymptotes = max(0, n - m)

        # -------------------------------------------------------------
        # Rule 1: Number of Branches & Terminations
        # -------------------------------------------------------------
        rule1_summary = (
            f"Rule 1: Open-loop transfer function has n = {n} poles and m = {m} zeros. "
            f"The root locus has {num_branches} branches; {num_asymptotes} branches terminate at infinity."
        )
        steps.append(rule1_summary)

        poles_fmt = ", ".join(
            f"{p.real:.4g}" if abs(p.imag) < 1e-6 else f"{p.real:.4g} ± {abs(p.imag):.4g}j"
            for p in parsed_poles
        )
        zeros_fmt = (
            ", ".join(
                f"{z.real:.4g}" if abs(z.imag) < 1e-6 else f"{z.real:.4g} ± {abs(z.imag):.4g}j"
                for z in parsed_zeros
            )
            if parsed_zeros
            else r"\text{None}"
        )

        detailed_derivations.append(
            rf"**Rule 1 (Branches & Terminations):** Open-loop system has $n = {n}$ poles (${poles_fmt}$) "
            rf"and $m = {m}$ zeros (${zeros_fmt}$). The root locus consists of **{num_branches} branches** starting at the poles ($K=0$). "
            rf"**{m} branches** terminate at finite zeros, and **{num_asymptotes} branches** terminate at infinity ($\infty$) as $K \to \infty$."
        )

        # -------------------------------------------------------------
        # Rule 2: Real-Axis Segments
        # -------------------------------------------------------------
        real_pts: list[float] = []
        for p in parsed_poles:
            if abs(p.imag) < 1e-6:
                real_pts.append(float(p.real))
        for z in parsed_zeros:
            if abs(z.imag) < 1e-6:
                real_pts.append(float(z.real))

        real_pts.sort(reverse=True)
        real_axis_segments: list[tuple[float, float]] = []

        for idx in range(len(real_pts)):
            if (idx + 1) % 2 == 1:
                high = real_pts[idx]
                low = real_pts[idx + 1] if idx + 1 < len(real_pts) else -float("inf")
                real_axis_segments.append((low, high))

        if real_axis_segments:
            seg_strs = []
            for low, high in real_axis_segments:
                low_str = r"-\infty" if math.isinf(low) else f"{low:.4g}"
                high_str = r"\infty" if math.isinf(high) else f"{high:.4g}"
                seg_strs.append(rf"[{low_str}, {high_str}]")
            seg_tex = r" \cup ".join(seg_strs)
        else:
            seg_tex = r"\emptyset"

        rule2_summary = (
            f"Rule 2: Real-axis locus lies to the left of an odd number of real open-loop poles and zeros: "
            f"{real_axis_segments}."
        )
        steps.append(rule2_summary)
        detailed_derivations.append(
            rf"**Rule 2 (Real-Axis Segments):** A point on the real axis belongs to the root locus if and only if the total number of real open-loop poles and zeros to its right is odd: "
            rf"$$s \in {seg_tex}$$"
        )

        # -------------------------------------------------------------
        # Rule 3: Asymptotes (Centroid and Angles)
        # -------------------------------------------------------------
        centroid: float | None = None
        asymptote_angles: list[float] = []

        if num_asymptotes > 0:
            sum_poles = sum(p.real for p in parsed_poles)
            sum_zeros = sum(z.real for z in parsed_zeros)
            centroid = float((sum_poles - sum_zeros) / num_asymptotes)
            angles_formula_parts: list[str] = []
            for k in range(num_asymptotes):
                angle = ((2 * k + 1) * 180.0) / num_asymptotes
                norm_angle = round(_normalize_angle_deg(angle), 2)
                asymptote_angles.append(norm_angle)
                angles_formula_parts.append(
                    rf"\theta_{k} = \frac{{(2({k})+1) \cdot 180^\circ}}{{{num_asymptotes}}} = {norm_angle:g}^\circ"
                )

            rule3_summary = (
                f"Rule 3: Asymptote centroid sigma_a = ({sum_poles:.4g} - {sum_zeros:.4g}) / {num_asymptotes} = {centroid:.4g}. "
                f"Asymptote angles: {asymptote_angles}."
            )
            steps.append(rule3_summary)
            angles_joined = ", ".join(angles_formula_parts)
            detailed_derivations.append(
                rf"**Rule 3 (Asymptotes of Locus as $s \to \infty$):** "
                rf"The $n - m = {num_asymptotes}$ branches radiating toward infinity follow straight-line asymptotes intersecting the real axis at centroid $\sigma_a$:"
                rf"$$\sigma_a = \frac{{\sum_{{i=1}}^n \operatorname{{Re}}(p_i) - \sum_{{j=1}}^m \operatorname{{Re}}(z_j)}}{{n - m}} = \frac{{{sum_poles:.4g} - ({sum_zeros:.4g})}}{{{num_asymptotes}}} = {centroid:.4g}$$"
                rf"The asymptote angles with the positive real axis are:"
                rf"$$\theta_k = \frac{{(2k + 1) \cdot 180^\circ}}{{n - m}} \implies {angles_joined}$$"
            )
        else:
            steps.append("Rule 3: Number of asymptotes is 0 since n = m.")
            detailed_derivations.append(
                r"**Rule 3 (Asymptotes):** All branches terminate at finite open-loop zeros ($n = m \implies$ no asymptotes toward $\infty$)."
            )

        # -------------------------------------------------------------
        # Rule 4: Breakaway and Break-in Points via dK/ds = 0
        # -------------------------------------------------------------
        num_poly = np.poly1d(parsed_num)
        den_poly = np.poly1d(parsed_den)

        d_den = np.polyder(den_poly)
        d_num = np.polyder(num_poly)

        # Characteristic: D'(s)*N(s) - D(s)*N'(s) = 0
        p_ba = np.polysub(np.polymul(d_den, num_poly), np.polymul(den_poly, d_num))

        breakaway_points: list[dict[str, Any]] = []
        if len(p_ba) > 0 and not np.all(p_ba == 0):
            roots_ba = np.roots(p_ba)
            for r in roots_ba:
                if abs(r.imag) < 1e-5:
                    s_val = float(r.real)
                    num_val = float(np.polyval(parsed_num, s_val))
                    if abs(num_val) > 1e-12:
                        den_val = float(np.polyval(parsed_den, s_val))
                        k_val = -den_val / num_val
                        if k_val >= -1e-6:
                            k_val = max(0.0, k_val)
                            on_locus = any(
                                low - 1e-6 <= s_val <= high + 1e-6
                                for low, high in real_axis_segments
                            )
                            if on_locus:
                                delta = 1e-4
                                den_m = float(np.polyval(parsed_den, s_val - delta))
                                num_m = float(np.polyval(parsed_num, s_val - delta))
                                den_p = float(np.polyval(parsed_den, s_val + delta))
                                num_p = float(np.polyval(parsed_num, s_val + delta))

                                k_m = -den_m / num_m if abs(num_m) > 1e-12 else k_val
                                k_p = -den_p / num_p if abs(num_p) > 1e-12 else k_val

                                pt_type = (
                                    "breakaway"
                                    if (k_val >= k_m - 1e-6 and k_val >= k_p - 1e-6)
                                    else "break-in"
                                )

                                if not any(abs(bp["s"] - s_val) < 1e-4 for bp in breakaway_points):
                                    breakaway_points.append(
                                        {
                                            "s": round(s_val, 4),
                                            "k": round(k_val, 4),
                                            "type": pt_type,
                                        }
                                    )

        rule4_summary = (
            f"Rule 4: Solved dK/ds = 0 for breakaway/break-in candidates: {breakaway_points}."
        )
        steps.append(rule4_summary)

        if breakaway_points:
            bp_tex_list = [
                rf"s = {bp['s']:.4g} \text{{ ({bp['type']}, }} K = {bp['k']:.4g}\text{{)}}"
                for bp in breakaway_points
            ]
            detailed_derivations.append(
                rf"**Rule 4 (Breakaway / Break-in Points):** Solving $\frac{{dK}}{{ds}} = -\frac{{d}}{{ds}}\left[\frac{{D(s)}}{{N(s)}}\right] = 0 \implies D'(s)N(s) - D(s)N'(s) = 0$ yields valid real-axis locus points with $K \ge 0$:"
                rf"$$\text{{Breakaway/Break-in: }} {', '.join(bp_tex_list)}$$"
            )
        else:
            detailed_derivations.append(
                r"**Rule 4 (Breakaway / Break-in Points):** Solving $\frac{dK}{ds} = 0$ yields no valid points with $K \ge 0$ on the real-axis root locus segments."
            )

        # -------------------------------------------------------------
        # Rule 5: Departure Angles from Complex Poles
        # -------------------------------------------------------------
        departure_angles: dict[complex, float] = {}
        for idx, p in enumerate(parsed_poles):
            if abs(p.imag) > 1e-6:
                angle_p_sum = sum(
                    math.atan2(p.imag - other_p.imag, p.real - other_p.real)
                    for j, other_p in enumerate(parsed_poles)
                    if j != idx
                )
                angle_z_sum = sum(
                    math.atan2(p.imag - z.imag, p.real - z.real) for z in parsed_zeros
                )
                theta_dep_rad = math.pi - angle_p_sum + angle_z_sum
                theta_dep_deg = math.degrees(theta_dep_rad)
                departure_angles[p] = round(_normalize_angle_deg(theta_dep_deg), 2)

        if departure_angles:
            steps.append(
                f"Rule 5: Calculated departure angles for complex poles: {departure_angles}."
            )
            dep_tex_list = [
                rf"\theta_{{\text{{dep}}}}(p = {p.real:.3g} + {p.imag:.3g}j) = 180^\circ - \sum \theta_p + \sum \theta_z = {deg:.2f}^\circ"
                for p, deg in departure_angles.items()
                if p.imag > 0
            ]
            detailed_derivations.append(
                rf"**Rule 5 (Angles of Departure from Complex Poles):** "
                rf"By the phase condition $\angle G(s)H(s) = (2k+1)180^\circ$ at $s \approx p_i$:"
                rf"$$\theta_{{\text{{dep}}}} = 180^\circ - \sum \theta_{{p}} + \sum \theta_{{z}} \implies {', '.join(dep_tex_list)}$$"
            )
        else:
            detailed_derivations.append(
                r"**Rule 5 (Angles of Departure):** No complex open-loop poles exist ($\implies$ no departure angle calculations needed)."
            )

        # -------------------------------------------------------------
        # Rule 6: Arrival Angles to Complex Zeros
        # -------------------------------------------------------------
        arrival_angles: dict[complex, float] = {}
        for idx, z in enumerate(parsed_zeros):
            if abs(z.imag) > 1e-6:
                angle_z_sum = sum(
                    math.atan2(z.imag - other_z.imag, z.real - other_z.real)
                    for j, other_z in enumerate(parsed_zeros)
                    if j != idx
                )
                angle_p_sum = sum(
                    math.atan2(z.imag - p.imag, z.real - p.real) for p in parsed_poles
                )
                theta_arr_rad = math.pi - angle_z_sum + angle_p_sum
                theta_arr_deg = math.degrees(theta_arr_rad)
                arrival_angles[z] = round(_normalize_angle_deg(theta_arr_deg), 2)

        if arrival_angles:
            steps.append(f"Rule 6: Calculated arrival angles for complex zeros: {arrival_angles}.")
            arr_tex_list = [
                rf"\theta_{{\text{{arr}}}}(z = {z.real:.3g} + {z.imag:.3g}j) = 180^\circ - \sum \theta_z + \sum \theta_p = {deg:.2f}^\circ"
                for z, deg in arrival_angles.items()
                if z.imag > 0
            ]
            detailed_derivations.append(
                rf"**Rule 6 (Angles of Arrival at Complex Zeros):** "
                rf"By the phase condition $\angle G(s)H(s) = (2k+1)180^\circ$ at $s \approx z_j$:"
                rf"$$\theta_{{\text{{arr}}}} = 180^\circ - \sum \theta_{{z}} + \sum \theta_{{p}} \implies {', '.join(arr_tex_list)}$$"
            )
        else:
            detailed_derivations.append(
                r"**Rule 6 (Angles of Arrival):** No complex open-loop zeros exist ($\implies$ no arrival angle calculations needed)."
            )

        # -------------------------------------------------------------
        # Rule 7: Imaginary Axis Crossings (s = j*omega)
        # -------------------------------------------------------------
        imag_crossings: list[dict[str, Any]] = []
        w_scan = np.logspace(-3, 4, 3000)
        jw_scan = 1j * w_scan

        d_vals = np.polyval(parsed_den, jw_scan)
        n_vals = np.polyval(parsed_num, jw_scan)

        k_complex = -d_vals / n_vals
        k_imag = np.imag(k_complex)

        zero_crossings = np.where(np.diff(np.sign(k_imag)))[0]
        for idx_val in zero_crossings:
            idx_int = int(idx_val)
            w1, w2 = float(w_scan[idx_int]), float(w_scan[idx_int + 1])
            im1, im2 = float(k_imag[idx_int]), float(k_imag[idx_int + 1])
            if im2 != im1:
                w_root = float(w1 - im1 * (w2 - w1) / (im2 - im1))
                for _ in range(8):
                    d_val_c = complex(np.polyval(parsed_den, 1j * w_root))
                    n_val_c = complex(np.polyval(parsed_num, 1j * w_root))
                    if abs(n_val_c) < 1e-12:
                        break
                    k_val_c = -d_val_c / n_val_c
                    im_val = float(np.imag(k_val_c))
                    if abs(im_val) < 1e-9:
                        break
                    dw = 1e-6
                    k_val_dw = -complex(np.polyval(parsed_den, 1j * (w_root + dw))) / complex(
                        np.polyval(parsed_num, 1j * (w_root + dw))
                    )
                    d_im = float((np.imag(k_val_dw) - im_val) / dw)
                    if abs(d_im) < 1e-12:
                        break
                    w_root = float(w_root - im_val / d_im)

                d_cross = complex(np.polyval(parsed_den, 1j * w_root))
                n_cross = complex(np.polyval(parsed_num, 1j * w_root))
                k_at_cross = float(np.real(-d_cross / n_cross))

                if (
                    k_at_cross > 1e-6
                    and w_root > 1e-4
                    and not any(abs(c["omega"] - w_root) < 1e-3 for c in imag_crossings)
                ):
                    imag_crossings.append(
                        {
                            "omega": round(float(w_root), 4),
                            "k": round(k_at_cross, 4),
                            "s": complex(0, round(float(w_root), 4)),
                        }
                    )

        rule7_summary = f"Rule 7: Imaginary axis crossings found: {imag_crossings}."
        steps.append(rule7_summary)

        if imag_crossings:
            cross_tex_list = [
                rf"s = \pm j{c['omega']:.4g} \text{{ (}}\omega = {c['omega']:.4g}\text{{ rad/s) at }} K_{{\text{{crit}}}} = {c['k']:.4g}"
                for c in imag_crossings
            ]
            detailed_derivations.append(
                rf"**Rule 7 (Imaginary Axis Crossings & Marginal Stability):** "
                rf"Substituting $s = j\omega$ into the characteristic equation $D(j\omega) + K N(j\omega) = 0$ yields critical stability crossings with the $j\omega$-axis:"
                rf"$$\text{{Crossings: }} {', '.join(cross_tex_list)}$$"
            )
        else:
            detailed_derivations.append(
                r"**Rule 7 (Imaginary Axis Crossings):** The root locus branches do not intersect the imaginary $j\omega$-axis for any real $K > 0$."
            )

        self.num_poles = n
        self.num_zeros = m
        self.num_branches = num_branches
        self.num_asymptotes = num_asymptotes
        self.poles = parsed_poles
        self.zeros = parsed_zeros
        self.real_axis_segments = real_axis_segments
        self.centroid = centroid
        self.asymptote_angles_deg = asymptote_angles
        self.breakaway_points = breakaway_points
        self.departure_angles = departure_angles
        self.arrival_angles = arrival_angles
        self.imag_axis_crossings = imag_crossings
        self.steps = steps
        self._detailed_derivations = detailed_derivations

    def explain_steps(self) -> list[str]:
        """Return a structured step-by-step list of Markdown/LaTeX explanations for all Evans rules.

        Returns
        -------
        list[str]
            List of step-by-step markdown explanations detailing each rule derivation.
        """
        return list(self._detailed_derivations)

    def _repr_latex_(self) -> str:
        """Render the complete Root Locus derivations as formatted LaTeX for Jupyter."""
        lines = [r"\begin{aligned}"]
        lines.append(r"\textbf{Analytical Root Locus Derivation Summary:}" + r" \\")

        # Rule 1
        lines.append(
            rf"\textbf{{Rule 1 (Branches):}} &\quad n = {self.num_poles}\text{{ poles}}, "
            rf"m = {self.num_zeros}\text{{ zeros}} \implies {self.num_branches}\text{{ branches}}, "
            rf"{self.num_asymptotes}\text{{ asymptotes to }}\infty \\"
        )

        # Rule 2: Real axis
        if self.real_axis_segments:
            seg_strs = []
            for low, high in self.real_axis_segments:
                low_str = r"-\infty" if math.isinf(low) else f"{low:.3g}"
                high_str = r"\infty" if math.isinf(high) else f"{high:.3g}"
                seg_strs.append(rf"[{low_str}, {high_str}]")
            seg_tex = r" \cup ".join(seg_strs)
        else:
            seg_tex = r"\emptyset"
        lines.append(rf"\textbf{{Rule 2 (Real-Axis Segments):}} &\quad s \in {seg_tex} \\")

        # Rule 3: Asymptotes
        if self.num_asymptotes > 0 and self.centroid is not None:
            angles_str = ", ".join(rf"{deg:g}^\circ" for deg in self.asymptote_angles_deg)
            lines.append(
                rf"\textbf{{Rule 3 (Asymptotes):}} &\quad \sigma_a = {self.centroid:.4g}, \quad \theta_k \in \{{{angles_str}\}} \\"
            )
        else:
            lines.append(
                r"\textbf{Rule 3 (Asymptotes):} &\quad \text{No asymptotes to } \infty (n = m) \\"
            )

        # Rule 4: Breakaway points
        if self.breakaway_points:
            bp_strs = [
                rf"s = {bp['s']:.4g} \text{{ ({bp['type']}, }} K = {bp['k']:.4g}\text{{)}}"
                for bp in self.breakaway_points
            ]
            lines.append(rf"\textbf{{Rule 4 (Breakaway/Break-in):}} &\quad {', '.join(bp_strs)} \\")
        else:
            lines.append(
                r"\textbf{Rule 4 (Breakaway/Break-in):} &\quad \text{None on the root locus} \\"
            )

        # Rule 5: Departure angles
        if self.departure_angles:
            dep_strs = [
                rf"\theta_{{\text{{dep}}}}(p={p.real:.3g}{p.imag:+.3g}j) = {deg:.2f}^\circ"
                for p, deg in self.departure_angles.items()
                if p.imag > 0
            ]
            lines.append(rf"\textbf{{Rule 5 (Departure Angles):}} &\quad {', '.join(dep_strs)} \\")

        # Rule 6: Arrival angles
        if self.arrival_angles:
            arr_strs = [
                rf"\theta_{{\text{{arr}}}}(z={z.real:.3g}{z.imag:+.3g}j) = {deg:.2f}^\circ"
                for z, deg in self.arrival_angles.items()
                if z.imag > 0
            ]
            lines.append(rf"\textbf{{Rule 6 (Arrival Angles):}} &\quad {', '.join(arr_strs)} \\")

        # Rule 7: Imaginary axis crossings
        if self.imag_axis_crossings:
            cross_strs = [
                rf"s = \pm j{c['omega']:.4g} \text{{ at }} K_{{\text{{crit}}}} = {c['k']:.4g}"
                for c in self.imag_axis_crossings
            ]
            lines.append(
                rf"\textbf{{Rule 7 (j\omega Crossings):}} &\quad {', '.join(cross_strs)} \\"
            )
        else:
            lines.append(
                r"\textbf{Rule 7 (j\omega Crossings):} &\quad \text{No crossings with the imaginary axis} \\"
            )

        lines.append(r"\end{aligned}")
        return "\n".join(lines)

    def _repr_markdown_(self) -> str:
        """Return a Markdown representation for Jupyter environments."""
        return f"$${self._repr_latex_()}$$"

    def __str__(self) -> str:
        """Format a human-readable ASCII representation of Root Locus rules."""
        lines = ["=== Analytical Root Locus Rules Summary ==="]
        lines.append(
            f"Poles (n={self.num_poles}): {[complex(round(p.real, 4), round(p.imag, 4)) for p in self.poles]}"
        )
        lines.append(
            f"Zeros (m={self.num_zeros}): {[complex(round(z.real, 4), round(z.imag, 4)) for z in self.zeros]}"
        )
        lines.append(
            f"Branches: {self.num_branches}, Asymptotes to Infinity: {self.num_asymptotes}"
        )
        lines.append("")

        # Real axis
        seg_strs = [f"[{low:g}, {high:g}]" for low, high in self.real_axis_segments]
        lines.append(f"Real-Axis Locus Segments: {' U '.join(seg_strs) if seg_strs else 'None'}")

        # Asymptotes
        if self.centroid is not None:
            angles_str = ", ".join(f"{deg:g}°" for deg in self.asymptote_angles_deg)
            lines.append(f"Asymptote Centroid: sigma_a = {self.centroid:.4f}")
            lines.append(f"Asymptote Angles: {angles_str}")

        # Breakaway
        if self.breakaway_points:
            lines.append("Breakaway / Break-in Points:")
            for bp in self.breakaway_points:
                lines.append(f"  - s = {bp['s']:.4f} ({bp['type']}) at K = {bp['k']:.4f}")
        else:
            lines.append("Breakaway / Break-in Points: None")

        # Departure / Arrival
        if self.departure_angles:
            lines.append("Departure Angles from Complex Poles:")
            for p, deg in self.departure_angles.items():
                lines.append(f"  - Pole {p}: {deg:.2f}°")
        if self.arrival_angles:
            lines.append("Arrival Angles to Complex Zeros:")
            for z, deg in self.arrival_angles.items():
                lines.append(f"  - Zero {z}: {deg:.2f}°")

        # jw Crossings
        if self.imag_axis_crossings:
            lines.append("Imaginary Axis Crossings:")
            for c in self.imag_axis_crossings:
                lines.append(f"  - s = ±j{c['omega']:.4f} at K_crit = {c['k']:.4f}")
        else:
            lines.append("Imaginary Axis Crossings: None")

        if self.steps:
            lines.append("\nStep-by-Step Derivation:")
            for s_note in self.steps:
                lines.append(f"  * {s_note}")

        return "\n".join(lines)

explain_steps

explain_steps()

Return a structured step-by-step list of Markdown/LaTeX explanations for all Evans rules.

Returns:

Type	Description
list[str]	List of step-by-step markdown explanations detailing each rule derivation.

Source code in src/ctrlpy/symbolic/root_locus.py

def explain_steps(self) -> list[str]:
    """Return a structured step-by-step list of Markdown/LaTeX explanations for all Evans rules.

    Returns
    -------
    list[str]
        List of step-by-step markdown explanations detailing each rule derivation.
    """
    return list(self._detailed_derivations)

__str__

__str__()

Format a human-readable ASCII representation of Root Locus rules.

Source code in src/ctrlpy/symbolic/root_locus.py

def __str__(self) -> str:
    """Format a human-readable ASCII representation of Root Locus rules."""
    lines = ["=== Analytical Root Locus Rules Summary ==="]
    lines.append(
        f"Poles (n={self.num_poles}): {[complex(round(p.real, 4), round(p.imag, 4)) for p in self.poles]}"
    )
    lines.append(
        f"Zeros (m={self.num_zeros}): {[complex(round(z.real, 4), round(z.imag, 4)) for z in self.zeros]}"
    )
    lines.append(
        f"Branches: {self.num_branches}, Asymptotes to Infinity: {self.num_asymptotes}"
    )
    lines.append("")

    # Real axis
    seg_strs = [f"[{low:g}, {high:g}]" for low, high in self.real_axis_segments]
    lines.append(f"Real-Axis Locus Segments: {' U '.join(seg_strs) if seg_strs else 'None'}")

    # Asymptotes
    if self.centroid is not None:
        angles_str = ", ".join(f"{deg:g}°" for deg in self.asymptote_angles_deg)
        lines.append(f"Asymptote Centroid: sigma_a = {self.centroid:.4f}")
        lines.append(f"Asymptote Angles: {angles_str}")

    # Breakaway
    if self.breakaway_points:
        lines.append("Breakaway / Break-in Points:")
        for bp in self.breakaway_points:
            lines.append(f"  - s = {bp['s']:.4f} ({bp['type']}) at K = {bp['k']:.4f}")
    else:
        lines.append("Breakaway / Break-in Points: None")

    # Departure / Arrival
    if self.departure_angles:
        lines.append("Departure Angles from Complex Poles:")
        for p, deg in self.departure_angles.items():
            lines.append(f"  - Pole {p}: {deg:.2f}°")
    if self.arrival_angles:
        lines.append("Arrival Angles to Complex Zeros:")
        for z, deg in self.arrival_angles.items():
            lines.append(f"  - Zero {z}: {deg:.2f}°")

    # jw Crossings
    if self.imag_axis_crossings:
        lines.append("Imaginary Axis Crossings:")
        for c in self.imag_axis_crossings:
            lines.append(f"  - s = ±j{c['omega']:.4f} at K_crit = {c['k']:.4f}")
    else:
        lines.append("Imaginary Axis Crossings: None")

    if self.steps:
        lines.append("\nStep-by-Step Derivation:")
        for s_note in self.steps:
            lines.append(f"  * {s_note}")

    return "\n".join(lines)

root_locus_rules

root_locus_rules(
    sys_or_poly=None,
    *,
    num=None,
    den=None,
    poles=None,
    zeros=None,
    gain=1.0,
)

Derive formal classroom Root Locus analytical rules step-by-step.

Parameters:

Name	Type	Description	Default
sys_or_poly	TransferFunction | Any	Open-loop transfer function \(G(s)H(s)\), LTI model, or symbolic expression.	None
num	Sequence[float] | None	Numerator polynomial coefficients.	None
den	Sequence[float] | None	Denominator polynomial coefficients.	None
poles	Sequence[complex | float] | None	Open-loop pole locations.	None
zeros	Sequence[complex | float] | None	Open-loop zero locations.	None
gain	float	Open-loop scalar gain factor.	1.0

Returns:

Type	Description
RootLocusRules	Analyzed rules object with metrics, derivation steps, and LaTeX formatting.

Source code in src/ctrlpy/symbolic/root_locus.py

def root_locus_rules(
    sys_or_poly: Any = None,
    *,
    num: Sequence[float] | None = None,
    den: Sequence[float] | None = None,
    poles: Sequence[complex | float] | None = None,
    zeros: Sequence[complex | float] | None = None,
    gain: float = 1.0,
) -> RootLocusRules:
    """Derive formal classroom Root Locus analytical rules step-by-step.

    Parameters
    ----------
    sys_or_poly : TransferFunction | Any, optional
        Open-loop transfer function $G(s)H(s)$, LTI model, or symbolic expression.
    num : Sequence[float] | None, default=None
        Numerator polynomial coefficients.
    den : Sequence[float] | None, default=None
        Denominator polynomial coefficients.
    poles : Sequence[complex | float] | None, default=None
        Open-loop pole locations.
    zeros : Sequence[complex | float] | None, default=None
        Open-loop zero locations.
    gain : float, default=1.0
        Open-loop scalar gain factor.

    Returns
    -------
    RootLocusRules
        Analyzed rules object with metrics, derivation steps, and LaTeX formatting.
    """
    return RootLocusRules(
        sys_or_poly=sys_or_poly,
        num=num,
        den=den,
        poles=poles,
        zeros=zeros,
        gain=gain,
    )

Steady-State Error & Static Constants

steady_state

Steady-state error analysis and static error constants.

SteadyStateResult dataclass

Container for steady-state error analysis and static constants.

Attributes:

Name	Type	Description
system_type	int	Number of open-loop integrators (poles at s = 0).
kp	float	Position error constant Kp = lim(s->0) G(s).
kv	float	Velocity error constant Kv = lim(s->0) s*G(s).
ka	float	Acceleration error constant Ka = lim(s->0) s^2*G(s).
ess_step	float	Steady-state tracking error for unit step input: 1 / (1 + Kp).
ess_ramp	float	Steady-state tracking error for unit ramp input: 1 / Kv.
ess_parabolic	float	Steady-state tracking error for unit parabolic input: 1 / Ka.
is_closed_loop_stable	bool	True if all unity feedback closed-loop poles lie strictly in the open LHP.
closed_loop_poles	list[complex]	The poles of the closed-loop transfer function T(s) = G(s) / (1 + G(s)).
steps	list[str]	Pedagogical step-by-step explanatory notes.

Source code in src/ctrlpy/pedagogy/steady_state.py

@dataclass
class SteadyStateResult:
    """Container for steady-state error analysis and static constants.

    Attributes
    ----------
    system_type : int
        Number of open-loop integrators (poles at s = 0).
    kp : float
        Position error constant Kp = lim(s->0) G(s).
    kv : float
        Velocity error constant Kv = lim(s->0) s*G(s).
    ka : float
        Acceleration error constant Ka = lim(s->0) s^2*G(s).
    ess_step : float
        Steady-state tracking error for unit step input: 1 / (1 + Kp).
    ess_ramp : float
        Steady-state tracking error for unit ramp input: 1 / Kv.
    ess_parabolic : float
        Steady-state tracking error for unit parabolic input: 1 / Ka.
    is_closed_loop_stable : bool
        True if all unity feedback closed-loop poles lie strictly in the open LHP.
    closed_loop_poles : list[complex]
        The poles of the closed-loop transfer function T(s) = G(s) / (1 + G(s)).
    steps : list[str]
        Pedagogical step-by-step explanatory notes.
    """

    system_type: int
    kp: float
    kv: float
    ka: float
    ess_step: float
    ess_ramp: float
    ess_parabolic: float
    is_closed_loop_stable: bool
    closed_loop_poles: list[complex] = field(default_factory=list)
    steps: list[str] = field(default_factory=list)

    def _repr_latex_(self) -> str:
        """Render steady-state error metrics and table as formatted LaTeX for Jupyter."""

        def _fmt(val: float) -> str:
            if math.isinf(val):
                return r"\infty"
            if abs(val) < 1e-12:
                return "0"
            return f"{val:.4g}"

        cl_status = (
            r"\text{Stable (Final Value Theorem valid)}"
            if self.is_closed_loop_stable
            else r"\textbf{\textcolor{red}{Unstable (FVT invalid - steady-state errors diverge!)}}"
        )

        latex_str = (
            r"\begin{aligned}"
            rf"\textbf{{System Classification:}} &\quad \text{{Type {self.system_type}}} \\"
            rf"\textbf{{Closed-Loop Stability:}} &\quad {cl_status} \\"
            r"\end{aligned}" + "\n\n"
            r"\begin{array}{|l|c|c|}"
            r"\hline"
            r"\textbf{Test Input } r(t) & \textbf{Static Error Constant} & \textbf{Steady-State Error } e_{ss} \\"
            r"\hline"
            rf"\text{{Unit Step: }} 1(t) & K_p = {_fmt(self.kp)} & e_{{ss,\text{{step}}}} = {_fmt(self.ess_step)} \\"
            rf"\text{{Unit Ramp: }} t \cdot 1(t) & K_v = {_fmt(self.kv)} & e_{{ss,\text{{ramp}}}} = {_fmt(self.ess_ramp)} \\"
            rf"\text{{Unit Parabola: }} \frac{{1}}{{2}}t^2 \cdot 1(t) & K_a = {_fmt(self.ka)} & e_{{ss,\text{{parabolic}}}} = {_fmt(self.ess_parabolic)} \\"
            r"\hline"
            r"\end{array}"
        )
        return latex_str

    def _repr_markdown_(self) -> str:
        """Return a Markdown representation for Jupyter environments."""
        return f"$${self._repr_latex_()}$$"

    def __str__(self) -> str:
        """Format a human-readable ASCII representation of steady-state analysis."""

        def _fmt(val: float) -> str:
            if math.isinf(val):
                return "inf"
            if abs(val) < 1e-12:
                return "0.0"
            return f"{val:.4f}"

        lines = ["=== Steady-State Error Analysis ==="]
        lines.append(f"System Type: Type {self.system_type}")
        lines.append(f"Closed-Loop Asymptotically Stable: {self.is_closed_loop_stable}")
        if not self.is_closed_loop_stable:
            lines.append(
                "  [WARNING: Closed-loop is unstable! Final Value Theorem does not apply.]"
            )
        lines.append("")
        lines.append(
            f"{'Input Type':<20} | {'Static Constant':<18} | {'Steady-State Error ess':<22}"
        )
        lines.append("-" * 66)
        lines.append(
            f"{'Unit Step 1(t)':<20} | {f'Kp = {_fmt(self.kp)}':<18} | {_fmt(self.ess_step):<22}"
        )
        lines.append(
            f"{'Unit Ramp t':<20} | {f'Kv = {_fmt(self.kv)}':<18} | {_fmt(self.ess_ramp):<22}"
        )
        lines.append(
            f"{'Unit Parabola 0.5*t^2':<20} | {f'Ka = {_fmt(self.ka)}':<18} | {_fmt(self.ess_parabolic):<22}"
        )
        lines.append("-" * 66)
        if self.steps:
            lines.append("\nAnalysis Notes:")
            for s in self.steps:
                lines.append(f"  * {s}")
        return "\n".join(lines)

__str__

__str__()

Format a human-readable ASCII representation of steady-state analysis.

Source code in src/ctrlpy/pedagogy/steady_state.py

def __str__(self) -> str:
    """Format a human-readable ASCII representation of steady-state analysis."""

    def _fmt(val: float) -> str:
        if math.isinf(val):
            return "inf"
        if abs(val) < 1e-12:
            return "0.0"
        return f"{val:.4f}"

    lines = ["=== Steady-State Error Analysis ==="]
    lines.append(f"System Type: Type {self.system_type}")
    lines.append(f"Closed-Loop Asymptotically Stable: {self.is_closed_loop_stable}")
    if not self.is_closed_loop_stable:
        lines.append(
            "  [WARNING: Closed-loop is unstable! Final Value Theorem does not apply.]"
        )
    lines.append("")
    lines.append(
        f"{'Input Type':<20} | {'Static Constant':<18} | {'Steady-State Error ess':<22}"
    )
    lines.append("-" * 66)
    lines.append(
        f"{'Unit Step 1(t)':<20} | {f'Kp = {_fmt(self.kp)}':<18} | {_fmt(self.ess_step):<22}"
    )
    lines.append(
        f"{'Unit Ramp t':<20} | {f'Kv = {_fmt(self.kv)}':<18} | {_fmt(self.ess_ramp):<22}"
    )
    lines.append(
        f"{'Unit Parabola 0.5*t^2':<20} | {f'Ka = {_fmt(self.ka)}':<18} | {_fmt(self.ess_parabolic):<22}"
    )
    lines.append("-" * 66)
    if self.steps:
        lines.append("\nAnalysis Notes:")
        for s in self.steps:
            lines.append(f"  * {s}")
    return "\n".join(lines)

steady_state_analysis

steady_state_analysis(sys)

Perform steady-state tracking error analysis for canonical test inputs.

Identifies system type (number of poles at s = 0), computes static error constants (Kp, Kv, Ka), evaluates steady-state errors (ess for step, ramp, parabola), and verifies closed-loop stability to validate the Final Value Theorem.

Parameters:

Name	Type	Description	Default
sys	TransferFunction	The open-loop LTI transfer function G(s) (assuming unity negative feedback).	required

Returns:

Type	Description
SteadyStateResult	Dataclass containing the system type, static error constants, steady-state errors, and closed-loop stability verification.

Examples:

>>> import ctrlpy as cp
>>> from ctrlpy.pedagogy import steady_state_analysis
>>> G = cp.tf([10], [1, 2, 0])  # Type 1 system
>>> res = steady_state_analysis(G)
>>> res.system_type
1
>>> res.kp
inf
>>> res.kv
5.0
>>> res.ess_step
0.0
>>> res.ess_ramp
0.2

Source code in src/ctrlpy/pedagogy/steady_state.py

def steady_state_analysis(sys: TransferFunction | Any) -> SteadyStateResult:
    """Perform steady-state tracking error analysis for canonical test inputs.

    Identifies system type (number of poles at s = 0), computes static error constants
    (Kp, Kv, Ka), evaluates steady-state errors (ess for step, ramp, parabola), and
    verifies closed-loop stability to validate the Final Value Theorem.

    Parameters
    ----------
    sys : TransferFunction
        The open-loop LTI transfer function G(s) (assuming unity negative feedback).

    Returns
    -------
    SteadyStateResult
        Dataclass containing the system type, static error constants, steady-state errors,
        and closed-loop stability verification.

    Examples
    --------
    >>> import ctrlpy as cp
    >>> from ctrlpy.pedagogy import steady_state_analysis
    >>> G = cp.tf([10], [1, 2, 0])  # Type 1 system
    >>> res = steady_state_analysis(G)
    >>> res.system_type
    1
    >>> res.kp
    inf
    >>> res.kv
    5.0
    >>> res.ess_step
    0.0
    >>> res.ess_ramp
    0.2
    """
    if not isinstance(sys, TransferFunction):
        if hasattr(sys, "to_tf"):
            sys = sys.to_tf()
        else:
            raise TypeError(
                "Expected a TransferFunction instance or an LTI object convertible to TransferFunction."
            )

    steps: list[str] = []

    # 1. Identify System Type (number of poles at the origin s = 0)
    den = list(sys.den)
    # Count trailing zeros in denominator
    type_n = 0
    for coeff in reversed(den):
        if abs(coeff) < 1e-12:
            type_n += 1
        else:
            break

    steps.append(
        f"Step 1: Open-loop denominator has {type_n} trailing zero(s) at s = 0 -> Classified as Type {type_n} system."
    )

    # 2. Compute Static Error Constants
    # G(s) = N(s) / (s^N * D_bar(s))
    # D_bar(0) is the last non-zero coefficient of den
    d_bar_0 = den[len(den) - 1 - type_n] if (len(den) - 1 - type_n) >= 0 else 1.0
    n_0 = sys.num[-1]

    k_bode = float(n_0 / d_bar_0)

    # Position Error Constant: Kp = lim_{s->0} G(s)
    if type_n == 0:
        kp = k_bode
    else:
        kp = float("inf")

    # Velocity Error Constant: Kv = lim_{s->0} s*G(s)
    if type_n == 0:
        kv = 0.0
    elif type_n == 1:
        kv = k_bode
    else:
        kv = float("inf")

    # Acceleration Error Constant: Ka = lim_{s->0} s^2*G(s)
    if type_n < 2:
        ka = 0.0
    elif type_n == 2:
        ka = k_bode
    else:
        ka = float("inf")

    steps.append(
        f"Step 2: Static error constants: Kp = {kp:g}, Kv = {kv:g}, Ka = {ka:g} (Bode gain K_0 = {k_bode:.4g})."
    )

    # 3. Compute Steady-State Errors ess
    # Step input
    if math.isinf(kp):
        ess_step = 0.0
    else:
        ess_step = 1.0 / (1.0 + kp)

    # Ramp input
    if kv == 0.0:
        ess_ramp = float("inf")
    elif math.isinf(kv):
        ess_ramp = 0.0
    else:
        ess_ramp = 1.0 / kv

    # Parabolic input
    if ka == 0.0:
        ess_parabolic = float("inf")
    elif math.isinf(ka):
        ess_parabolic = 0.0
    else:
        ess_parabolic = 1.0 / ka

    steps.append(
        f"Step 3: Steady-state errors: ess(step) = {ess_step:g}, ess(ramp) = {ess_ramp:g}, ess(parabola) = {ess_parabolic:g}."
    )

    # 4. Verify Closed-Loop Stability
    # Characteristic equation D_cl(s) = D(s) + N(s)
    pad_len = max(len(sys.den), len(sys.num))
    den_padded = np.pad(sys.den, (pad_len - len(sys.den), 0))
    num_padded = np.pad(sys.num, (pad_len - len(sys.num), 0))
    cl_den = den_padded + num_padded

    cl_roots = np.roots(cl_den)
    cl_poles = [complex(r) for r in cl_roots]
    is_cl_stable = all(p.real < -1e-7 for p in cl_poles)

    if is_cl_stable:
        steps.append(
            "Step 4: Unity feedback closed-loop system is asymptotically stable (all closed-loop poles in LHP)."
        )
    else:
        steps.append(
            "Step 4: WARNING: Closed-loop system has unstable or imaginary poles! Final Value Theorem calculations are invalid."
        )

    return SteadyStateResult(
        system_type=type_n,
        kp=kp,
        kv=kv,
        ka=ka,
        ess_step=ess_step,
        ess_ramp=ess_ramp,
        ess_parabolic=ess_parabolic,
        is_closed_loop_stable=is_cl_stable,
        closed_loop_poles=cl_poles,
        steps=steps,
    )

State-Space Canonical Forms & Controllability/Observability Tutor

state_space

Analytical State-Space Canonical Forms and Controllability/Observability Tutor.

StateSpaceTutor

Pedagogical state-space canonical forms and controllability/observability tutor.

Provides rigorous analytical, step-by-step transformations and checks for LTI state-space systems: 1. Controllability Analysis: Symbolically builds \(\mathcal{C} = [B \quad AB \dots A^{n-1}B]\), computes determinant, rank, and nullspace. 2. Observability Analysis: Symbolically builds \(\mathcal{O} = [C^T \quad A^TC^T \dots (A^{n-1})^TC^T]^T\), computes determinant, rank, and unobservable subspaces. 3. Popov-Belevitch-Hautus (PBH) Mode Decomposition: Evaluates eigenvalue-by-eigenvalue rank tests \(\operatorname{rank}[\lambda_i I - A \quad B]\) and \(\operatorname{rank}[\lambda_i I - A; C]\) to classify all modes into the 4 Kalman canonical subspaces (\(\Sigma_{co}, \Sigma_{c\bar{o}}, \Sigma_{\bar{c}o}, \Sigma_{\bar{c}\bar{o}}\)). 4. Controllable Canonical Form (CCF / Phase-Variable): Derives \((A_c, B_c, C_c, D)\) and transformation \(T_c = \mathcal{C}\mathcal{W}\). 5. Observable Canonical Form (OCF): Derives \((A_o, B_o, C_o, D)\) and transformation \(T_o = \mathcal{O}^{-1}\mathcal{C}_c\). 6. Jordan / Diagonal Modal Canonical Form: Derives modal matrix \(V\) (eigenvectors/generalized eigenvectors) and decoupled state equations. 7. Step-by-Step Pedagogical Derivations: .explain_steps() returns an exhaustive mathematical trace.

Parameters:

Name	Type	Description	Default
sys_or_A	Any	StateSpace model, TransferFunction model, \(A\) matrix, or polynomial numerator.	None
B	Any	Input matrix \(B \in \mathbb{R}^{n \times m}\).	None
C	Any	Output matrix \(C \in \mathbb{R}^{p \times n}\).	None
D	Any	Direct feedthrough matrix \(D \in \mathbb{R}^{p \times m}\).	None
num	Sequence[Any] | None	Transfer function numerator polynomial coefficients if initializing from transfer function.	None
den	Sequence[Any] | None	Transfer function denominator polynomial coefficients if initializing from transfer function.	None
dt	float | None	Sampling period for discrete-time systems.	None

Attributes:

Name	Type	Description
A	Matrix	State transition matrix (\(n \times n\)).
B	Matrix	Input matrix (\(n \times m\)).
C	Matrix	Output matrix (\(p \times n\)).
D	Matrix	Direct transmission matrix (\(p \times m\)).
n_states	int	Number of state variables \(n\).
n_inputs	int	Number of inputs \(m\).
n_outputs	int	Number of outputs \(p\).
dt	float | None	Sampling period (None for continuous-time).
controllability_matrix	Matrix	The Kalman controllability matrix \(\mathcal{C}\).
observability_matrix	Matrix	The Kalman observability matrix \(\mathcal{O}\).
controllability_rank	int	Rank of \(\mathcal{C}\).
observability_rank	int	Rank of \(\mathcal{O}\).
is_controllable	bool	Whether the system is fully controllable (\(\operatorname{rank}(\mathcal{C}) = n\)).
is_observable	bool	Whether the system is fully observable (\(\operatorname{rank}(\mathcal{O}) = n\)).
characteristic_polynomial	Poly	Symbolic characteristic polynomial \(p(s) = \det(sI - A)\).
eigenvalues	list[Expr]	List of all eigenvalues \(\lambda_i\).
modes	list[ModeAnalysis]	Detailed PBH and Kalman modal classification for each eigenvalue.
uncontrollable_modes	list[Expr]	List of eigenvalues that fail PBH controllability.
unobservable_modes	list[Expr]	List of eigenvalues that fail PBH observability.
transfer_function	Expr	Analytical transfer function \(G(s) = C(sI - A)^{-1}B + D\).
ccf	CanonicalFormResult	Controllable Canonical Form transformation result.
ocf	CanonicalFormResult	Observable Canonical Form transformation result.
jcf	CanonicalFormResult	Jordan / Diagonal Canonical Form transformation result.
steps	list[str]	Pedagogical step-by-step mathematical explanations.

Source code in src/ctrlpy/symbolic/state_space.py

class StateSpaceTutor:
    r"""Pedagogical state-space canonical forms and controllability/observability tutor.

    Provides rigorous analytical, step-by-step transformations and checks for LTI state-space systems:
    1. **Controllability Analysis**: Symbolically builds $\mathcal{C} = [B \quad AB \dots A^{n-1}B]$,
       computes determinant, rank, and nullspace.
    2. **Observability Analysis**: Symbolically builds $\mathcal{O} = [C^T \quad A^TC^T \dots (A^{n-1})^TC^T]^T$,
       computes determinant, rank, and unobservable subspaces.
    3. **Popov-Belevitch-Hautus (PBH) Mode Decomposition**: Evaluates eigenvalue-by-eigenvalue rank tests
       $\operatorname{rank}[\lambda_i I - A \quad B]$ and $\operatorname{rank}[\lambda_i I - A; C]$ to classify
       all modes into the 4 Kalman canonical subspaces ($\Sigma_{co}, \Sigma_{c\bar{o}}, \Sigma_{\bar{c}o}, \Sigma_{\bar{c}\bar{o}}$).
    4. **Controllable Canonical Form (CCF / Phase-Variable)**: Derives $(A_c, B_c, C_c, D)$ and transformation $T_c = \mathcal{C}\mathcal{W}$.
    5. **Observable Canonical Form (OCF)**: Derives $(A_o, B_o, C_o, D)$ and transformation $T_o = \mathcal{O}^{-1}\mathcal{C}_c$.
    6. **Jordan / Diagonal Modal Canonical Form**: Derives modal matrix $V$ (eigenvectors/generalized eigenvectors) and decoupled state equations.
    7. **Step-by-Step Pedagogical Derivations**: `.explain_steps()` returns an exhaustive mathematical trace.

    Parameters
    ----------
    sys_or_A : Any, optional
        StateSpace model, TransferFunction model, $A$ matrix, or polynomial numerator.
    B : Any, optional
        Input matrix $B \in \mathbb{R}^{n \times m}$.
    C : Any, optional
        Output matrix $C \in \mathbb{R}^{p \times n}$.
    D : Any, optional
        Direct feedthrough matrix $D \in \mathbb{R}^{p \times m}$.
    num : Sequence[Any] | None, default=None
        Transfer function numerator polynomial coefficients if initializing from transfer function.
    den : Sequence[Any] | None, default=None
        Transfer function denominator polynomial coefficients if initializing from transfer function.
    dt : float | None, default=None
        Sampling period for discrete-time systems.

    Attributes
    ----------
    A : sp.Matrix
        State transition matrix ($n \times n$).
    B : sp.Matrix
        Input matrix ($n \times m$).
    C : sp.Matrix
        Output matrix ($p \times n$).
    D : sp.Matrix
        Direct transmission matrix ($p \times m$).
    n_states : int
        Number of state variables $n$.
    n_inputs : int
        Number of inputs $m$.
    n_outputs : int
        Number of outputs $p$.
    dt : float | None
        Sampling period (None for continuous-time).
    controllability_matrix : sp.Matrix
        The Kalman controllability matrix $\mathcal{C}$.
    observability_matrix : sp.Matrix
        The Kalman observability matrix $\mathcal{O}$.
    controllability_rank : int
        Rank of $\mathcal{C}$.
    observability_rank : int
        Rank of $\mathcal{O}$.
    is_controllable : bool
        Whether the system is fully controllable ($\operatorname{rank}(\mathcal{C}) = n$).
    is_observable : bool
        Whether the system is fully observable ($\operatorname{rank}(\mathcal{O}) = n$).
    characteristic_polynomial : sp.Poly
        Symbolic characteristic polynomial $p(s) = \det(sI - A)$.
    eigenvalues : list[sp.Expr]
        List of all eigenvalues $\lambda_i$.
    modes : list[ModeAnalysis]
        Detailed PBH and Kalman modal classification for each eigenvalue.
    uncontrollable_modes : list[sp.Expr]
        List of eigenvalues that fail PBH controllability.
    unobservable_modes : list[sp.Expr]
        List of eigenvalues that fail PBH observability.
    transfer_function : sp.Expr
        Analytical transfer function $G(s) = C(sI - A)^{-1}B + D$.
    ccf : CanonicalFormResult
        Controllable Canonical Form transformation result.
    ocf : CanonicalFormResult
        Observable Canonical Form transformation result.
    jcf : CanonicalFormResult
        Jordan / Diagonal Canonical Form transformation result.
    steps : list[str]
        Pedagogical step-by-step mathematical explanations.
    """

    def __init__(
        self,
        sys_or_A: Any = None,
        B: Any = None,
        C: Any = None,
        D: Any = None,
        *,
        num: Sequence[Any] | None = None,
        den: Sequence[Any] | None = None,
        dt: float | None = None,
    ) -> None:
        self.dt = dt
        steps: list[str] = []

        # 1. Parse Inputs into SymPy Matrices
        A_sp, B_sp, C_sp, D_sp, inferred_dt = self._parse_inputs(sys_or_A, B, C, D, num, den)
        if self.dt is None:
            self.dt = inferred_dt

        self.A = A_sp
        self.B = B_sp
        self.C = C_sp
        self.D = D_sp

        n = self.A.rows
        m = self.B.cols
        p = self.C.rows

        self.n_states = n
        self.n_inputs = m
        self.n_outputs = p

        sym_var = sp.Symbol("s") if self.dt is None else sp.Symbol("z")
        var_str = "s" if self.dt is None else "z"

        # 2. Compute Characteristic Polynomial and Eigenvalues
        sI_minus_A = sym_var * sp.eye(n) - self.A
        char_poly_expr = sp.simplify(sI_minus_A.det())
        char_poly = sp.Poly(char_poly_expr, sym_var)
        self.char_poly_expr = char_poly_expr
        self.characteristic_polynomial = char_poly

        # Monic polynomial coefficients
        coeffs = char_poly.all_coeffs()
        lead = coeffs[0]
        monic_coeffs = [sp.simplify(c / lead) for c in coeffs]
        # p(s) = s^n + a_{n-1} s^{n-1} + ... + a_1 s + a_0
        # monic_coeffs has length n+1: [1, a_{n-1}, a_{n-2}, ..., a_0]
        a_coeffs = monic_coeffs[1:]  # [a_{n-1}, a_{n-2}, ..., a_0]
        # a_dict: a_dict[0] = a_0, a_dict[1] = a_1, ..., a_dict[n-1] = a_{n-1}
        a_dict: dict[int, sp.Expr] = {}
        for idx, coef in enumerate(reversed(a_coeffs)):
            a_dict[idx] = coef

        eigen_dict = self.A.eigenvals()
        eigenvalues: list[sp.Expr] = []
        eigenvalue_multiplicities: dict[sp.Expr, int] = {}
        for ev, mult in eigen_dict.items():
            eigenvalues.append(ev)
            eigenvalue_multiplicities[ev] = mult
        self.eigenvalues = eigenvalues
        self.eigenvalue_multiplicities = eigenvalue_multiplicities

        # 3. Compute Controllability Matrix
        ctrl_cols: list[sp.Matrix] = []
        intermediate_ab: list[sp.Matrix] = []
        for k in range(n):
            a_pow = self.A**k
            ab_k = a_pow * self.B
            intermediate_ab.append(ab_k)
            ctrl_cols.append(ab_k)

        self.controllability_matrix = sp.Matrix.hstack(*ctrl_cols) if ctrl_cols else sp.zeros(n, 0)
        self.controllability_rank = int(self.controllability_matrix.rank())
        self.is_controllable = bool(self.controllability_rank == n)

        # 4. Compute Observability Matrix
        obs_rows: list[sp.Matrix] = []
        intermediate_ca: list[sp.Matrix] = []
        for k in range(n):
            a_pow = self.A**k
            ca_k = self.C * a_pow
            intermediate_ca.append(ca_k)
            obs_rows.append(ca_k)

        self.observability_matrix = sp.Matrix.vstack(*obs_rows) if obs_rows else sp.zeros(0, n)
        self.observability_rank = int(self.observability_matrix.rank())
        self.is_observable = bool(self.observability_rank == n)

        # 5. Compute Analytical Transfer Function
        try:
            adj_matrix = sI_minus_A.adjugate()
            tf_matrix = sp.simplify(
                (self.C * adj_matrix * self.B + self.D * char_poly_expr) / char_poly_expr
            )
            if self.is_siso:
                self.transfer_function = sp.simplify(tf_matrix[0, 0])
                num_expr, den_expr = sp.fraction(self.transfer_function)
                self.tf_numerator = sp.Poly(num_expr, sym_var)
                self.tf_denominator = sp.Poly(den_expr, sym_var)
            else:
                self.transfer_function = tf_matrix
                self.tf_numerator = None
                self.tf_denominator = None
        except (TypeError, ValueError, AttributeError, sp.SympifyError, sp.MatrixError):
            self.transfer_function = sp.sympify(0)
            self.tf_numerator = None
            self.tf_denominator = None

        # 6. PBH Eigenvalue Mode Analysis
        modes: list[ModeAnalysis] = []
        uncontrollable_modes: list[sp.Expr] = []
        unobservable_modes: list[sp.Expr] = []

        for ev in eigenvalues:
            mult = eigenvalue_multiplicities[ev]
            # Controllability PBH: [ev * I - A, B]
            pbh_c_mat = sp.Matrix.hstack(ev * sp.eye(n) - self.A, self.B)
            c_rank = int(pbh_c_mat.rank())
            ev_ctrl = bool(c_rank == n)

            # Observability PBH: [ev * I - A; C]
            pbh_o_mat = sp.Matrix.vstack(ev * sp.eye(n) - self.A, self.C)
            o_rank = int(pbh_o_mat.rank())
            ev_obs = bool(o_rank == n)

            if not ev_ctrl:
                uncontrollable_modes.append(ev)
            if not ev_obs:
                unobservable_modes.append(ev)

            if ev_ctrl and ev_obs:
                k_type = "co"
                desc = f"Mode λ = {sp.latex(ev)} is both Controllable and Observable (present in transfer function)."
            elif ev_ctrl and not ev_obs:
                k_type = "c_unobs"
                desc = f"Mode λ = {sp.latex(ev)} is Controllable but Unobservable (actuated by input, hidden from output / pole-zero cancelled)."
            elif not ev_ctrl and ev_obs:
                k_type = "unctrl_obs"
                desc = f"Mode λ = {sp.latex(ev)} is Uncontrollable but Observable (visible in output, unaffected by input / pole-zero cancelled)."
            else:
                k_type = "unctrl_unobs"
                desc = f"Mode λ = {sp.latex(ev)} is Uncontrollable and Unobservable (completely decoupled internal dynamic)."

            modes.append(
                ModeAnalysis(
                    eigenvalue=ev,
                    multiplicity=mult,
                    is_controllable=ev_ctrl,
                    is_observable=ev_obs,
                    pbh_c_rank=c_rank,
                    pbh_o_rank=o_rank,
                    kalman_type=k_type,
                    description=desc,
                )
            )

        self.modes = modes
        self.uncontrollable_modes = uncontrollable_modes
        self.unobservable_modes = unobservable_modes

        # 7. Canonical Form Transformations (SISO Companion Forms & Jordan Modal Form)
        self.ccf = self._derive_ccf(a_dict, sym_var)
        self.ocf = self._derive_ocf(a_dict, sym_var)
        self.jcf = self._derive_jcf()

        # 8. Compile Step-by-Step Pedagogical Derivations
        self._build_pedagogical_steps(steps, var_str, a_dict, intermediate_ab, intermediate_ca)
        self.steps = steps

    @property
    def is_siso(self) -> bool:
        """Return True if system is Single-Input Single-Output."""
        return bool(self.n_inputs == 1 and self.n_outputs == 1)

    @staticmethod
    def _parse_inputs(
        sys_or_A: Any,
        B: Any,
        C: Any,
        D: Any,
        num: Sequence[Any] | None,
        den: Sequence[Any] | None,
    ) -> tuple[sp.Matrix, sp.Matrix, sp.Matrix, sp.Matrix, float | None]:
        """Parse various system representations into standard SymPy (A, B, C, D) matrices."""
        inferred_dt: float | None = None

        # Case 1: Transfer function numerator / denominator sequences provided
        if num is not None and den is not None:
            return StateSpaceTutor._from_tf_coeffs(num, den, dt=inferred_dt)

        # Case 2: TransferFunction model instance
        if hasattr(sys_or_A, "num") and hasattr(sys_or_A, "den"):
            inferred_dt = getattr(sys_or_A, "dt", None)
            return StateSpaceTutor._from_tf_coeffs(sys_or_A.num, sys_or_A.den, dt=inferred_dt)

        # Case 3: StateSpace model instance
        if (
            hasattr(sys_or_A, "A")
            and hasattr(sys_or_A, "B")
            and hasattr(sys_or_A, "C")
            and hasattr(sys_or_A, "D")
        ):
            inferred_dt = getattr(sys_or_A, "dt", None)
            return (
                _to_sympy_matrix(sys_or_A.A),
                _to_sympy_matrix(sys_or_A.B),
                _to_sympy_matrix(sys_or_A.C),
                _to_sympy_matrix(sys_or_A.D),
                inferred_dt,
            )

        # Case 4: Explicit matrices A, B, C, D
        if sys_or_A is not None:
            a_sp = _to_sympy_matrix(sys_or_A)
            if a_sp.rows != a_sp.cols:
                raise DimensionMismatchError(
                    f"Matrix A must be square, got shape ({a_sp.rows}, {a_sp.cols})."
                )
            n = a_sp.rows

            if B is None:
                raise ValueError("Matrix B must be provided when matrix A is given.")
            b_sp = _to_sympy_matrix(B)
            if b_sp.rows != n:
                if b_sp.cols == n and b_sp.rows == 1:
                    b_sp = b_sp.T
                elif b_sp.rows == 1 and n == 1:
                    pass
                else:
                    raise DimensionMismatchError(
                        f"Matrix B rows ({b_sp.rows}) must match matrix A dimension ({n})."
                    )

            if C is None:
                raise ValueError("Matrix C must be provided when matrix A is given.")
            c_sp = _to_sympy_matrix(C)
            if c_sp.cols != n:
                if c_sp.rows == n and c_sp.cols == 1:
                    c_sp = c_sp.T
                elif c_sp.cols == 1 and n == 1:
                    pass
                else:
                    raise DimensionMismatchError(
                        f"Matrix C columns ({c_sp.cols}) must match matrix A dimension ({n})."
                    )

            m = b_sp.cols
            p = c_sp.rows

            if D is None:
                d_sp = sp.zeros(p, m)
            else:
                d_sp = _to_sympy_matrix(D)
                if d_sp.rows != p or d_sp.cols != m:
                    if d_sp.rows == 1 and d_sp.cols == 1 and (p != 1 or m != 1):
                        d_sp = sp.Matrix(
                            [[d_sp[0, 0] if i == j else 0 for j in range(m)] for i in range(p)]
                        )
                    else:
                        raise DimensionMismatchError(
                            f"Matrix D shape ({d_sp.rows}, {d_sp.cols}) must match outputs {p} and inputs {m}."
                        )

            return a_sp, b_sp, c_sp, d_sp, inferred_dt

        raise ValueError(
            "Must provide either a StateSpace/TransferFunction model or explicit A, B, C, D matrices."
        )

    @staticmethod
    def _from_tf_coeffs(
        num: Sequence[Any],
        den: Sequence[Any],
        dt: float | None = None,
    ) -> tuple[sp.Matrix, sp.Matrix, sp.Matrix, sp.Matrix, float | None]:
        """Derive standard controllable companion state-space matrices from transfer function coefficients."""
        s = sp.Symbol("s")
        clean_num = [_to_sympy_matrix(x)[0, 0] for x in num]
        clean_den = [_to_sympy_matrix(x)[0, 0] for x in den]

        # Build polynomials
        num_poly = sp.Poly(clean_num, s)
        den_poly = sp.Poly(clean_den, s)

        den_coeffs = den_poly.all_coeffs()
        lead = den_coeffs[0]
        monic_den_coeffs = [sp.simplify(c / lead) for c in den_coeffs]
        monic_num_coeffs = [sp.simplify(c / lead) for c in num_poly.all_coeffs()]

        n = len(monic_den_coeffs) - 1
        if n <= 0:
            # Scalar gain
            d_val = monic_num_coeffs[0] if monic_num_coeffs else sp.Integer(0)
            return sp.Matrix([[0]]), sp.Matrix([[0]]), sp.Matrix([[0]]), sp.Matrix([[d_val]]), dt

        # Pad numerator to length n + 1
        num_padded = [sp.Integer(0)] * (n + 1 - len(monic_num_coeffs)) + monic_num_coeffs
        d_val = num_padded[0]

        # Strictly proper numerator: N'(s) = N(s) - D * Den(s)
        beta_coeffs = [
            sp.simplify(num_padded[i] - d_val * monic_den_coeffs[i]) for i in range(1, n + 1)
        ]
        # beta_coeffs: [beta_{n-1}, beta_{n-2}, ..., beta_0]

        # CCF Companion matrix A:
        # [ 0, 1, 0, ..., 0 ]
        # [ 0, 0, 1, ..., 0 ]
        # [ -a_0, -a_1, ..., -a_{n-1} ]
        a_den = list(reversed(monic_den_coeffs[1:]))  # [a_0, a_1, ..., a_{n-1}]
        a_rows: list[list[sp.Expr]] = []
        for i in range(n - 1):
            row = [sp.Integer(0)] * n
            row[i + 1] = sp.Integer(1)
            a_rows.append(row)
        a_rows.append([-a_den[i] for i in range(n)])

        A_sp = _to_sympy_matrix(sp.Matrix(a_rows))
        B_sp = sp.zeros(n, 1)
        B_sp[n - 1, 0] = sp.Integer(1)

        # C = [beta_0, beta_1, ..., beta_{n-1}]
        c_entries = list(reversed(beta_coeffs))  # [beta_0, beta_1, ..., beta_{n-1}]
        C_sp = _to_sympy_matrix(sp.Matrix([c_entries]))
        D_sp = _to_sympy_matrix(sp.Matrix([[d_val]]))

        return A_sp, B_sp, C_sp, D_sp, dt

    def _derive_ccf(self, a_dict: dict[int, sp.Expr], sym_var: sp.Symbol) -> CanonicalFormResult:
        """Derive Controllable Canonical Form (Phase-Variable Form) and transformation matrix."""
        n = self.n_states
        if not self.is_siso:
            return CanonicalFormResult(
                self.A,
                self.B,
                self.C,
                self.D,
                None,
                None,
                form_type="controllable",
                is_valid=False,
                explanation="Controllable Canonical Form (Phase-Variable) is defined for SISO systems.",
                dt=self.dt,
            )

        # Build companion matrices (Ac, Bc)
        # Ac has 1s on superdiagonal and bottom row [-a_0, -a_1, ..., -a_{n-1}]
        a_rows: list[list[sp.Expr]] = []
        for i in range(n - 1):
            row = [sp.Integer(0)] * n
            row[i + 1] = sp.Integer(1)
            a_rows.append(row)
        a_rows.append([-a_dict[i] for i in range(n)])
        Ac = sp.Matrix(a_rows)

        Bc = sp.zeros(n, 1)
        Bc[n - 1, 0] = sp.Integer(1)

        # Compute strictly proper transfer function numerator coefficients beta_0 ... beta_{n-1}
        # G(s) = C(sI - A)^-1 B + D
        tf_strictly_proper = sp.simplify(self.transfer_function - self.D[0, 0])
        # Multiply strictly proper transfer function by characteristic polynomial to obtain full uncancelled numerator
        num_uncancelled = sp.simplify(tf_strictly_proper * self.char_poly_expr)
        num_poly = sp.Poly(num_uncancelled, sym_var)
        num_coeffs = num_poly.all_coeffs()
        padded_num = [sp.Integer(0)] * (n - len(num_coeffs)) + num_coeffs
        # padded_num is [beta_{n-1}, beta_{n-2}, ..., beta_0]
        c_entries = list(reversed(padded_num))  # [beta_0, beta_1, ..., beta_{n-1}]
        Cc = sp.Matrix([c_entries])
        Dc = self.D

        if not self.is_controllable:
            return CanonicalFormResult(
                Ac,
                Bc,
                Cc,
                Dc,
                None,
                None,
                form_type="controllable",
                is_valid=False,
                explanation=(
                    f"System is uncontrollable (rank(C) = {self.controllability_rank} < {n}). "
                    "A nonsingular similarity transformation to Controllable Canonical Form does not exist."
                ),
                dt=self.dt,
            )

        # Compute controllability matrix of companion form: Cc_mat = [Bc, Ac*Bc, ..., Ac^(n-1)*Bc]
        cols_cc: list[sp.Matrix] = []
        for k in range(n):
            cols_cc.append((Ac**k) * Bc)
        Cc_mat = sp.Matrix.hstack(*cols_cc)

        # Transformation matrix: Tc = C * Cc_mat^-1 = C * W
        Cc_inv = Cc_mat.inv()
        Tc = sp.simplify(self.controllability_matrix * Cc_inv)
        Tc_inv = sp.simplify(Tc.inv())

        return CanonicalFormResult(
            Ac,
            Bc,
            Cc,
            Dc,
            Tc,
            Tc_inv,
            form_type="controllable",
            is_valid=True,
            explanation=(
                "Controllable Canonical Form (Phase-Variable form) derived via similarity transformation "
                r"$x = T_c z_c$, where $T_c = \mathcal{C} \mathcal{C}_c^{-1} = \mathcal{C} \mathcal{W}$."
            ),
            dt=self.dt,
        )

    def _derive_ocf(self, a_dict: dict[int, sp.Expr], sym_var: sp.Symbol) -> CanonicalFormResult:
        """Derive Observable Canonical Form and transformation matrix."""
        n = self.n_states
        if not self.is_siso:
            return CanonicalFormResult(
                self.A,
                self.B,
                self.C,
                self.D,
                None,
                None,
                form_type="observable",
                is_valid=False,
                explanation="Observable Canonical Form is defined for SISO systems.",
                dt=self.dt,
            )

        # Companion Ao = Ac^T
        # Rows:
        # [0, 0, ..., -a_0]
        # [1, 0, ..., -a_1]
        # [0, 1, ..., -a_2]
        # [0, 0, ..., -a_{n-1}]
        ao_rows: list[list[sp.Expr]] = []
        for i in range(n):
            row = [sp.Integer(0)] * n
            if i > 0:
                row[i - 1] = sp.Integer(1)
            row[n - 1] = -a_dict[i]
            ao_rows.append(row)
        Ao = sp.Matrix(ao_rows)

        # Strictly proper numerator coefficients
        tf_strictly_proper = sp.simplify(self.transfer_function - self.D[0, 0])
        num_uncancelled = sp.simplify(tf_strictly_proper * self.char_poly_expr)
        num_poly = sp.Poly(num_uncancelled, sym_var)
        num_coeffs = num_poly.all_coeffs()
        padded_num = [sp.Integer(0)] * (n - len(num_coeffs)) + num_coeffs
        # [beta_{n-1}, ..., beta_0] -> Bo = [beta_0, beta_1, ..., beta_{n-1}]^T
        b_entries = list(reversed(padded_num))
        Bo = sp.Matrix([[b] for b in b_entries])

        # Co = [0, 0, ..., 1]
        Co = sp.zeros(1, n)
        Co[0, n - 1] = sp.Integer(1)
        Do = self.D

        if not self.is_observable:
            return CanonicalFormResult(
                Ao,
                Bo,
                Co,
                Do,
                None,
                None,
                form_type="observable",
                is_valid=False,
                explanation=(
                    f"System is unobservable (rank(O) = {self.observability_rank} < {n}). "
                    "A nonsingular similarity transformation to Observable Canonical Form does not exist."
                ),
                dt=self.dt,
            )

        # Observability matrix of companion form: Oo_mat = [Co; Co*Ao; ...; Co*Ao^(n-1)]
        rows_oo: list[sp.Matrix] = []
        for k in range(n):
            rows_oo.append(Co * (Ao**k))
        Oo_mat = sp.Matrix.vstack(*rows_oo)

        # Transformation matrix: O * To = Oo_mat => To = O^-1 * Oo_mat
        O_inv = self.observability_matrix.inv()
        To = sp.simplify(O_inv * Oo_mat)
        To_inv = sp.simplify(To.inv())

        return CanonicalFormResult(
            Ao,
            Bo,
            Co,
            Do,
            To,
            To_inv,
            form_type="observable",
            is_valid=True,
            explanation=(
                "Observable Canonical Form derived via similarity transformation "
                r"$x = T_o z_o$, where $T_o = \mathcal{O}^{-1} \mathcal{O}_o$."
            ),
            dt=self.dt,
        )

    def _derive_jcf(self) -> CanonicalFormResult:
        """Derive Diagonal / Jordan Canonical Form and modal transformation matrix."""
        try:
            P, J = self.A.jordan_form()
            P = sp.simplify(P)
            J = sp.simplify(J)
            P_inv = sp.simplify(P.inv())

            Ad = J
            Bd = sp.simplify(P_inv * self.B)
            Cd = sp.simplify(self.C * P)
            Dd = self.D

            return CanonicalFormResult(
                Ad,
                Bd,
                Cd,
                Dd,
                P,
                P_inv,
                form_type="jordan",
                is_valid=True,
                explanation=(
                    "Jordan / Diagonal Canonical Form derived via modal transformation "
                    r"$x = V z_d$, where columns of $V$ are eigenvectors and generalized eigenvectors."
                ),
                dt=self.dt,
            )
        except (TypeError, ValueError, AttributeError, sp.SympifyError, sp.MatrixError) as err:
            return CanonicalFormResult(
                self.A,
                self.B,
                self.C,
                self.D,
                None,
                None,
                form_type="jordan",
                is_valid=False,
                explanation=f"Failed to compute Jordan canonical form: {err}",
                dt=self.dt,
            )

    def _build_pedagogical_steps(
        self,
        steps: list[str],
        var_str: str,
        a_dict: dict[int, sp.Expr],
        intermediate_ab: list[sp.Matrix],
        intermediate_ca: list[sp.Matrix],
    ) -> None:
        """Assemble structured LaTeX/Markdown explanations for all derivation stages."""
        n = self.n_states
        # 1. System Representation & Characteristic Polynomial
        poly_str = " + ".join(
            f"{sp.latex(a_dict[i])}{var_str}^{{{i}}}"
            if i > 1
            else (f"{sp.latex(a_dict[i])}{var_str}" if i == 1 else f"{sp.latex(a_dict[i])}")
            for i in reversed(range(n))
            if a_dict[i] != 0
        )
        if poly_str:
            poly_full = f"{var_str}^{{{n}}} + {poly_str}"
        else:
            poly_full = f"{var_str}^{{{n}}}"

        ev_items = [
            f"{sp.latex(ev)} \\text{{ (mult: {mult})}}"
            for ev, mult in self.eigenvalue_multiplicities.items()
        ]
        ev_str = ", ".join(ev_items)

        steps.append(
            f"### Step 1: System Definition & Characteristic Polynomial\n"
            f"The continuous/discrete state-space matrices are defined as:\n\n"
            f"$$A = {_sp_matrix_to_latex(self.A)}, \\quad B = {_sp_matrix_to_latex(self.B)}, "
            f"\\quad C = {_sp_matrix_to_latex(self.C)}, \\quad D = {_sp_matrix_to_latex(self.D)}$$\n\n"
            f"The characteristic polynomial is obtained by evaluating $p({var_str}) = \\det({var_str} I - A)$:\n\n"
            f"$$p({var_str}) = \\det\\left({var_str} I - A\\right) = {sp.latex(self.char_poly_expr)} = {poly_full} = 0$$\n\n"
            f"The open-loop eigenvalues (system poles) are:\n\n"
            f"$$\\lambda = \\left\\{{ {ev_str} \\right\\}}$$"
        )

        # 2. Controllability Derivation
        ab_terms = " \\quad ".join(
            rf"A^{k} B = {_sp_matrix_to_latex(mat)}" for k, mat in enumerate(intermediate_ab)
        )
        ctrl_det_str = (
            rf"\det(\mathcal{{C}}) = {sp.latex(self.controllability_matrix.det())}"
            if self.controllability_matrix.rows == self.controllability_matrix.cols
            else r"\text{N/A (Non-square matrix)}"
        )
        if self.is_controllable:
            ctrl_verdict = (
                r"The system is **fully controllable** ($\operatorname{rank}(\mathcal{C}) = n$). "
                r"All state trajectories can be steered arbitrarily in finite time."
            )
        else:
            ctrl_verdict = (
                rf"The system is **uncontrollable** ($\operatorname{{rank}}(\mathcal{{C}}) = {self.controllability_rank} < {n}$). "
                rf"A {n - self.controllability_rank}-dimensional uncontrollable subspace exists."
            )

        steps.append(
            f"### Step 2: Controllability Matrix & Kalman Rank Test\n"
            f"The Kalman controllability matrix $\\mathcal{{C}} \\in \\mathbb{{R}}^{{{n} \\times {self.controllability_matrix.cols}}}$ "
            f"is constructed from successive matrix powers $A^k B$:\n\n"
            f"$$\\mathcal{{C}} = \\begin{{bmatrix}} B & AB & A^2 B & \\dots & A^{{{n - 1}}} B \\end{{bmatrix}}$$\n\n"
            f"Evaluating individual block columns:\n\n"
            f"$${ab_terms}$$\n\n"
            f"Concatenating yields the full controllability matrix:\n\n"
            f"$$\\mathcal{{C}} = {_sp_matrix_to_latex(self.controllability_matrix)}$$\n\n"
            f"Determinant and rank evaluation:\n\n"
            f"$${ctrl_det_str}, \\quad \\operatorname{{rank}}(\\mathcal{{C}}) = {self.controllability_rank} \\quad (n = {n})$$\n\n"
            f"**Controllability Verdict:** {ctrl_verdict}"
        )

        # 3. Observability Derivation
        ca_terms = " \\quad ".join(
            rf"C A^{k} = {_sp_matrix_to_latex(mat)}" for k, mat in enumerate(intermediate_ca)
        )
        obs_det_str = (
            rf"\det(\mathcal{{O}}) = {sp.latex(self.observability_matrix.det())}"
            if self.observability_matrix.rows == self.observability_matrix.cols
            else r"\text{N/A (Non-square matrix)}"
        )

        if self.is_observable:
            obs_verdict = (
                r"The system is **fully observable** ($\operatorname{rank}(\mathcal{O}) = n$). "
                r"The initial state $x(0)$ can be uniquely reconstructed from output history $y(t)$."
            )
        else:
            obs_verdict = (
                rf"The system is **unobservable** ($\operatorname{{rank}}(\mathcal{{O}}) = {self.observability_rank} < {n}$). "
                rf"A {n - self.observability_rank}-dimensional unobservable subspace exists."
            )

        steps.append(
            f"### Step 3: Observability Matrix & Kalman Rank Test\n"
            f"The Kalman observability matrix $\\mathcal{{O}} \\in \\mathbb{{R}}^{{{self.observability_matrix.rows} \\times {n}}}$ "
            f"is constructed from successive output matrix products $C A^k$:\n\n"
            f"$$\\mathcal{{O}} = \\begin{{bmatrix}} C \\\\ CA \\\\ CA^2 \\\\ \\vdots \\\\ CA^{{{n - 1}}} \\end{{bmatrix}}$$\n\n"
            f"Evaluating individual block rows:\n\n"
            f"$${ca_terms}$$\n\n"
            f"Concatenating yields the full observability matrix:\n\n"
            f"$$\\mathcal{{O}} = {_sp_matrix_to_latex(self.observability_matrix)}$$\n\n"
            f"Determinant and rank evaluation:\n\n"
            f"$${obs_det_str}, \\quad \\operatorname{{rank}}(\\mathcal{{O}}) = {self.observability_rank} \\quad (n = {n})$$\n\n"
            f"**Observability Verdict:** {obs_verdict}"
        )

        # 4. PBH Eigenvalue Test & Mode Decomposition
        mode_rows: list[str] = []
        for m_obj in self.modes:
            mode_rows.append(
                f"- **Mode $\\lambda = {sp.latex(m_obj.eigenvalue)}$** (Algebraic Mult: {m_obj.multiplicity}):\n"
                f"  - Controllability PBH: $\\operatorname{{rank}}\\begin{{bmatrix}} \\lambda I - A & B \\end{{bmatrix}} = {m_obj.pbh_c_rank} / {n} "
                f"\\implies \\textbf{{{('Controllable' if m_obj.is_controllable else 'Uncontrollable')}}}$\n"
                f"  - Observability PBH: $\\operatorname{{rank}}\\begin{{bmatrix}} \\lambda I - A \\\\ C \\end{{bmatrix}} = {m_obj.pbh_o_rank} / {n} "
                f"\\implies \\textbf{{{('Observable' if m_obj.is_observable else 'Unobservable')}}}$\n"
                f"  - Kalman Subspace Classification: `{m_obj.kalman_type.upper()}` — {m_obj.description}"
            )
        mode_text = "\n".join(mode_rows)

        steps.append(
            f"### Step 4: Popov-Belevitch-Hautus (PBH) Mode Decomposition\n"
            f"The PBH rank tests check the modal controllability and observability of each eigenvalue $\\lambda_i$:\n\n"
            f"$$\\operatorname{{rank}}\\begin{{bmatrix}} \\lambda_i I - A & B \\end{{bmatrix}} = n \\iff \\text{{Mode }} \\lambda_i \\text{{ Controllable}}$$\n"
            f"$$\\operatorname{{rank}}\\begin{{bmatrix}} \\lambda_i I - A \\\\ C \\end{{bmatrix}} = n \\iff \\text{{Mode }} \\lambda_i \\text{{ Observable}}$$\n\n"
            f"{mode_text}"
        )

        # 5. Controllable Canonical Form
        if self.ccf.is_valid and self.ccf.T is not None:
            steps.append(
                f"### Step 5: Controllable Canonical Form (Phase-Variable)\n"
                f"The phase-variable Controllable Canonical Form $(A_c, B_c, C_c, D)$ is given by:\n\n"
                f"$$A_c = {_sp_matrix_to_latex(self.ccf.A)}, \\quad B_c = {_sp_matrix_to_latex(self.ccf.B)}, "
                f"\\quad C_c = {_sp_matrix_to_latex(self.ccf.C)}, \\quad D = {_sp_matrix_to_latex(self.ccf.D)}$$\n\n"
                f"Transformation matrix $T_c$ ($x = T_c z_c$) and inverse $T_c^{{-1}}$ ($z_c = T_c^{{-1}} x$):\n\n"
                f"$$T_c = \\mathcal{{C}} \\mathcal{{C}}_c^{{-1}} = {_sp_matrix_to_latex(self.ccf.T)}, \\quad "
                f"T_c^{{-1}} = {_sp_matrix_to_latex(self.ccf.T_inv)}$$\n\n"
                f"Verification: $T_c^{{-1}} A T_c = A_c$, $T_c^{{-1}} B = B_c$, $C T_c = C_c$."
            )
        else:
            steps.append(
                f"### Step 5: Controllable Canonical Form (Phase-Variable)\n{self.ccf.explanation}"
            )

        # 6. Observable Canonical Form
        if self.ocf.is_valid and self.ocf.T is not None:
            steps.append(
                f"### Step 6: Observable Canonical Form\n"
                f"The Observable Canonical Form $(A_o, B_o, C_o, D)$ is dual to the controllable form:\n\n"
                f"$$A_o = {_sp_matrix_to_latex(self.ocf.A)}, \\quad B_o = {_sp_matrix_to_latex(self.ocf.B)}, "
                f"\\quad C_o = {_sp_matrix_to_latex(self.ocf.C)}, \\quad D = {_sp_matrix_to_latex(self.ocf.D)}$$\n\n"
                f"Transformation matrix $T_o$ ($x = T_o z_o$) and inverse $T_o^{{-1}}$ ($z_o = T_o^{{-1}} x$):\n\n"
                f"$$T_o = \\mathcal{{O}}^{{-1}} \\mathcal{{O}}_o = {_sp_matrix_to_latex(self.ocf.T)}, \\quad "
                f"T_o^{{-1}} = {_sp_matrix_to_latex(self.ocf.T_inv)}$$\n\n"
                f"Verification: $T_o^{{-1}} A T_o = A_o$, $T_o^{{-1}} B = B_o$, $C T_o = C_o$."
            )
        else:
            steps.append(f"### Step 6: Observable Canonical Form\n{self.ocf.explanation}")

        # 7. Jordan Canonical Form
        if self.jcf.is_valid and self.jcf.T is not None:
            steps.append(
                f"### Step 7: Jordan / Diagonal Modal Form\n"
                f"The modal transformation matrix $V$ diagonalizes the state matrix $A$ into Jordan form:\n\n"
                f"$$A_d = {_sp_matrix_to_latex(self.jcf.A)}, \\quad B_d = {_sp_matrix_to_latex(self.jcf.B)}, "
                f"\\quad C_d = {_sp_matrix_to_latex(self.jcf.C)}, \\quad D = {_sp_matrix_to_latex(self.jcf.D)}$$\n\n"
                f"Modal eigenvector matrix $V$ ($x = V z_d$) and inverse $V^{{-1}}$ ($z_d = V^{{-1}} x$):\n\n"
                f"$$V = {_sp_matrix_to_latex(self.jcf.T)}, \\quad V^{{-1}} = {_sp_matrix_to_latex(self.jcf.T_inv)}$$\n\n"
                f"Decoupled modal equations: $\\dot{{z}}_{{d,i}} = \\lambda_i z_{{d,i}} + b_{{d,i}} u$, $y = \\sum c_{{d,i}} z_{{d,i}} + D u$."
            )
        else:
            steps.append(f"### Step 7: Jordan / Diagonal Modal Form\n{self.jcf.explanation}")

        # 8. Transfer Function & Cancellations
        if self.is_siso:
            steps.append(
                f"### Step 8: Analytical Transfer Function & Pole-Zero Cancellation\n"
                f"The input-output transfer function $G({var_str}) = C({var_str} I - A)^{{-1}} B + D$ is:\n\n"
                f"$$G({var_str}) = {sp.latex(self.transfer_function)}$$\n\n"
                f"- Characteristic polynomial order: $n = {n}$\n"
                f"- Transfer function minimal denominator order: $\\deg(D({var_str})) = {self.tf_denominator.degree() if self.tf_denominator is not None else n}$\n"
                f"- Uncontrollable poles: {[sp.latex(p) for p in self.uncontrollable_modes] or 'None'}\n"
                f"- Unobservable poles: {[sp.latex(p) for p in self.unobservable_modes] or 'None'}\n"
                f"- Poles cancelled in transfer function: {f'Yes ({n - (self.tf_denominator.degree() if self.tf_denominator is not None else n)} cancelled mode(s))' if len(self.uncontrollable_modes) > 0 or len(self.unobservable_modes) > 0 else 'None (Minimal Realization)'}"
            )

    def explain_steps(self) -> list[str]:
        """Return the step-by-step list of Markdown/LaTeX explanations.

        Returns
        -------
        list[str]
            List of formatted derivation steps.
        """
        return list(self.steps)

    def controllable_canonical_form(self) -> CanonicalFormResult:
        """Return the Controllable Canonical Form (Phase-Variable) result.

        Returns
        -------
        CanonicalFormResult
            CCF state matrices, transformation matrix, and derivation status.
        """
        return self.ccf

    def observable_canonical_form(self) -> CanonicalFormResult:
        """Return the Observable Canonical Form result.

        Returns
        -------
        CanonicalFormResult
            OCF state matrices, transformation matrix, and derivation status.
        """
        return self.ocf

    def jordan_canonical_form(self) -> CanonicalFormResult:
        """Return the Jordan / Diagonal Canonical Form result.

        Returns
        -------
        CanonicalFormResult
            Jordan/diagonal modal state matrices, modal matrix V, and derivation status.
        """
        return self.jcf

    def pbh_test(self) -> list[dict[str, Any]]:
        """Return summary of PBH controllability and observability test for each eigenvalue.

        Returns
        -------
        list[dict[str, Any]]
            List of dictionaries containing eigenvalue details, PBH ranks, and Kalman types.
        """
        return [
            {
                "eigenvalue": m.eigenvalue,
                "multiplicity": m.multiplicity,
                "is_controllable": m.is_controllable,
                "is_observable": m.is_observable,
                "pbh_controllability_rank": m.pbh_c_rank,
                "pbh_observability_rank": m.pbh_o_rank,
                "kalman_type": m.kalman_type,
                "description": m.description,
            }
            for m in self.modes
        ]

    def to_ss(self) -> StateSpace:
        """Convert system to numerical StateSpace model.

        Returns
        -------
        StateSpace
            Numerical StateSpace object.
        """
        from ctrlpy.models.state_space import StateSpace

        def _to_np(m: sp.Matrix) -> np.ndarray:
            arr = np.zeros((m.rows, m.cols), dtype=np.float64)
            for i in range(m.rows):
                for j in range(m.cols):
                    val = m[i, j]
                    if hasattr(val, "free_symbols") and val.free_symbols:
                        raise ValueError(
                            f"Cannot convert symbolic expression {val} to numeric float."
                        )
                    arr[i, j] = float(sp.re(val.evalf()))
            return arr

        return StateSpace(_to_np(self.A), _to_np(self.B), _to_np(self.C), _to_np(self.D))

    def to_tf(self) -> TransferFunction:
        """Convert system to numerical TransferFunction model.

        Returns
        -------
        TransferFunction
            Numerical TransferFunction object.

        Raises
        ------
        NotImplementedError
            If the system is MIMO.
        """
        from ctrlpy.models.transfer_function import TransferFunction

        if not self.is_siso:
            raise NotImplementedError("to_tf() is only supported for SISO systems.")

        if self.tf_numerator is None or self.tf_denominator is None:
            return self.to_ss().to_tf()

        num_coeffs = [float(sp.re(c.evalf())) for c in self.tf_numerator.all_coeffs()]
        den_coeffs = [float(sp.re(c.evalf())) for c in self.tf_denominator.all_coeffs()]
        return TransferFunction(num_coeffs, den_coeffs)

    def _repr_latex_(self) -> str:
        r"""Render a comprehensive LaTeX breakdown of the state-space system and canonical transformations."""
        dt_sym = r"\dot{x}" if self.dt is None else "x[k+1]"
        state_sym = "x" if self.dt is None else "x[k]"
        u_sym = "u" if self.dt is None else "u[k]"
        y_sym = "y" if self.dt is None else "y[k]"

        a_tex = _sp_matrix_to_latex(self.A)
        b_tex = _sp_matrix_to_latex(self.B)
        c_tex = _sp_matrix_to_latex(self.C)
        d_tex = _sp_matrix_to_latex(self.D)
        ctrl_tex = _sp_matrix_to_latex(self.controllability_matrix)
        obs_tex = _sp_matrix_to_latex(self.observability_matrix)

        ctrl_status = (
            r"\textbf{Controllable}" if self.is_controllable else r"\textbf{Uncontrollable}"
        )
        obs_status = r"\textbf{Observable}" if self.is_observable else r"\textbf{Unobservable}"

        ev_joined = ", ".join(sp.latex(ev) for ev in self.eigenvalues)
        lines = [
            r"\begin{aligned}",
            r"\textbf{State-Space Representation:} \\",
            rf"{dt_sym} &= {a_tex} {state_sym} + {b_tex} {u_sym} \\",
            rf"{y_sym} &= {c_tex} {state_sym} + {d_tex} {u_sym} \\[1em]",
            (
                r"\textbf{Controllability Matrix } (\mathcal{C}): \quad "
                rf"\operatorname{{rank}}(\mathcal{{C}}) = {self.controllability_rank} / {self.n_states} \implies {ctrl_status} \\"
            ),
            rf"\mathcal{{C}} &= {ctrl_tex} \\[1em]",
            (
                r"\textbf{Observability Matrix } (\mathcal{O}): \quad "
                rf"\operatorname{{rank}}(\mathcal{{O}}) = {self.observability_rank} / {self.n_states} \implies {obs_status} \\"
            ),
            rf"\mathcal{{O}} &= {obs_tex} \\[1em]",
            (
                r"\textbf{Characteristic Polynomial:} \quad "
                rf"p(s) = {sp.latex(self.char_poly_expr)} = 0 \\"
            ),
            (
                r"\textbf{Eigenvalues / Modes:} \quad "
                rf"\lambda = \left\{{ {ev_joined} \right\}}"
            ),
            r"\end{aligned}",
        ]
        return "\n".join(lines)

    def _repr_markdown_(self) -> str:
        """Render Markdown for Jupyter environments."""
        return f"$${self._repr_latex_()}$$"

    def __str__(self) -> str:
        """ASCII formatted breakdown of state-space system and canonical forms."""
        lines = [
            "=== Pedagogical State-Space Analysis ===",
            f"State Dimension: n={self.n_states}, Inputs: m={self.n_inputs}, Outputs: p={self.n_outputs}",
            f"A =\n{self.A}",
            f"B =\n{self.B}",
            f"C =\n{self.C}",
            f"D =\n{self.D}",
            "",
            f"Controllability Matrix C (Rank {self.controllability_rank}/{self.n_states}, Controllable={self.is_controllable}):",
            f"{self.controllability_matrix}",
            "",
            f"Observability Matrix O (Rank {self.observability_rank}/{self.n_states}, Observable={self.is_observable}):",
            f"{self.observability_matrix}",
            "",
            f"Characteristic Polynomial: {self.char_poly_expr} = 0",
            f"Eigenvalues: {self.eigenvalues}",
            f"Uncontrollable Modes: {self.uncontrollable_modes}",
            f"Unobservable Modes: {self.unobservable_modes}",
        ]
        return "\n".join(lines)

    def __repr__(self) -> str:
        return (
            f"StateSpaceTutor(n_states={self.n_states}, is_controllable={self.is_controllable}, "
            f"is_observable={self.is_observable})"
        )

is_siso property

is_siso

Return True if system is Single-Input Single-Output.

explain_steps

explain_steps()

Return the step-by-step list of Markdown/LaTeX explanations.

Returns:

Type	Description
list[str]	List of formatted derivation steps.

Source code in src/ctrlpy/symbolic/state_space.py

def explain_steps(self) -> list[str]:
    """Return the step-by-step list of Markdown/LaTeX explanations.

    Returns
    -------
    list[str]
        List of formatted derivation steps.
    """
    return list(self.steps)

controllable_canonical_form

controllable_canonical_form()

Return the Controllable Canonical Form (Phase-Variable) result.

Returns:

Type	Description
CanonicalFormResult	CCF state matrices, transformation matrix, and derivation status.

Source code in src/ctrlpy/symbolic/state_space.py

def controllable_canonical_form(self) -> CanonicalFormResult:
    """Return the Controllable Canonical Form (Phase-Variable) result.

    Returns
    -------
    CanonicalFormResult
        CCF state matrices, transformation matrix, and derivation status.
    """
    return self.ccf

observable_canonical_form

observable_canonical_form()

Return the Observable Canonical Form result.

Returns:

Type	Description
CanonicalFormResult	OCF state matrices, transformation matrix, and derivation status.

Source code in src/ctrlpy/symbolic/state_space.py

def observable_canonical_form(self) -> CanonicalFormResult:
    """Return the Observable Canonical Form result.

    Returns
    -------
    CanonicalFormResult
        OCF state matrices, transformation matrix, and derivation status.
    """
    return self.ocf

jordan_canonical_form

jordan_canonical_form()

Return the Jordan / Diagonal Canonical Form result.

Returns:

Type	Description
CanonicalFormResult	Jordan/diagonal modal state matrices, modal matrix V, and derivation status.

Source code in src/ctrlpy/symbolic/state_space.py

def jordan_canonical_form(self) -> CanonicalFormResult:
    """Return the Jordan / Diagonal Canonical Form result.

    Returns
    -------
    CanonicalFormResult
        Jordan/diagonal modal state matrices, modal matrix V, and derivation status.
    """
    return self.jcf

pbh_test

pbh_test()

Return summary of PBH controllability and observability test for each eigenvalue.

Returns:

Type	Description
list[dict[str, Any]]	List of dictionaries containing eigenvalue details, PBH ranks, and Kalman types.

Source code in src/ctrlpy/symbolic/state_space.py

def pbh_test(self) -> list[dict[str, Any]]:
    """Return summary of PBH controllability and observability test for each eigenvalue.

    Returns
    -------
    list[dict[str, Any]]
        List of dictionaries containing eigenvalue details, PBH ranks, and Kalman types.
    """
    return [
        {
            "eigenvalue": m.eigenvalue,
            "multiplicity": m.multiplicity,
            "is_controllable": m.is_controllable,
            "is_observable": m.is_observable,
            "pbh_controllability_rank": m.pbh_c_rank,
            "pbh_observability_rank": m.pbh_o_rank,
            "kalman_type": m.kalman_type,
            "description": m.description,
        }
        for m in self.modes
    ]

to_ss

to_ss()

Convert system to numerical StateSpace model.

Returns:

Type	Description
StateSpace	Numerical StateSpace object.

Source code in src/ctrlpy/symbolic/state_space.py

def to_ss(self) -> StateSpace:
    """Convert system to numerical StateSpace model.

    Returns
    -------
    StateSpace
        Numerical StateSpace object.
    """
    from ctrlpy.models.state_space import StateSpace

    def _to_np(m: sp.Matrix) -> np.ndarray:
        arr = np.zeros((m.rows, m.cols), dtype=np.float64)
        for i in range(m.rows):
            for j in range(m.cols):
                val = m[i, j]
                if hasattr(val, "free_symbols") and val.free_symbols:
                    raise ValueError(
                        f"Cannot convert symbolic expression {val} to numeric float."
                    )
                arr[i, j] = float(sp.re(val.evalf()))
        return arr

    return StateSpace(_to_np(self.A), _to_np(self.B), _to_np(self.C), _to_np(self.D))

to_tf

to_tf()

Convert system to numerical TransferFunction model.

Returns:

Type	Description
TransferFunction	Numerical TransferFunction object.

Raises:

Type	Description
NotImplementedError	If the system is MIMO.

Source code in src/ctrlpy/symbolic/state_space.py

def to_tf(self) -> TransferFunction:
    """Convert system to numerical TransferFunction model.

    Returns
    -------
    TransferFunction
        Numerical TransferFunction object.

    Raises
    ------
    NotImplementedError
        If the system is MIMO.
    """
    from ctrlpy.models.transfer_function import TransferFunction

    if not self.is_siso:
        raise NotImplementedError("to_tf() is only supported for SISO systems.")

    if self.tf_numerator is None or self.tf_denominator is None:
        return self.to_ss().to_tf()

    num_coeffs = [float(sp.re(c.evalf())) for c in self.tf_numerator.all_coeffs()]
    den_coeffs = [float(sp.re(c.evalf())) for c in self.tf_denominator.all_coeffs()]
    return TransferFunction(num_coeffs, den_coeffs)

__str__

__str__()

ASCII formatted breakdown of state-space system and canonical forms.

Source code in src/ctrlpy/symbolic/state_space.py

def __str__(self) -> str:
    """ASCII formatted breakdown of state-space system and canonical forms."""
    lines = [
        "=== Pedagogical State-Space Analysis ===",
        f"State Dimension: n={self.n_states}, Inputs: m={self.n_inputs}, Outputs: p={self.n_outputs}",
        f"A =\n{self.A}",
        f"B =\n{self.B}",
        f"C =\n{self.C}",
        f"D =\n{self.D}",
        "",
        f"Controllability Matrix C (Rank {self.controllability_rank}/{self.n_states}, Controllable={self.is_controllable}):",
        f"{self.controllability_matrix}",
        "",
        f"Observability Matrix O (Rank {self.observability_rank}/{self.n_states}, Observable={self.is_observable}):",
        f"{self.observability_matrix}",
        "",
        f"Characteristic Polynomial: {self.char_poly_expr} = 0",
        f"Eigenvalues: {self.eigenvalues}",
        f"Uncontrollable Modes: {self.uncontrollable_modes}",
        f"Unobservable Modes: {self.unobservable_modes}",
    ]
    return "\n".join(lines)

CanonicalFormResult

Analytical canonical form state-space system and transformation details.

Represents transformed state-space representations:

.. math::

\dot{z} = A_{\text{can}} z + B_{\text{can}} u, \quad y = C_{\text{can}} z + D u

under the state transformation \(x = T z\) (\(z = T^{-1} x\)).

Parameters:

Name	Type	Description	Default
A	Matrix	Canonical state matrix \(A_{\text{can}} = T^{-1} A T\).	required
B	Matrix	Canonical input matrix \(B_{\text{can}} = T^{-1} B\).	required
C	Matrix	Canonical output matrix \(C_{\text{can}} = C T\).	required
D	Matrix	Direct feedthrough matrix \(D\).	required
T	Matrix | None	Similarity transformation matrix \(T\) mapping canonical states \(z\) to original states \(x = T z\).	required
T_inv	Matrix | None	Inverse transformation matrix \(T^{-1}\) mapping original states to canonical states \(z = T^{-1} x\).	required
form_type	str	Type of canonical form ('controllable', 'observable', 'jordan').	required
is_valid	bool	Whether the canonical transformation is nonsingular and valid.	True
explanation	str	Analytical derivation and pedagogical notes for this canonical form.	''
dt	float | None	Sampling time for discrete systems (None for continuous-time).	None

Source code in src/ctrlpy/symbolic/state_space.py

class CanonicalFormResult:
    r"""Analytical canonical form state-space system and transformation details.

    Represents transformed state-space representations:

    .. math::

        \dot{z} = A_{\text{can}} z + B_{\text{can}} u, \quad y = C_{\text{can}} z + D u

    under the state transformation $x = T z$ ($z = T^{-1} x$).

    Parameters
    ----------
    A : sp.Matrix
        Canonical state matrix $A_{\text{can}} = T^{-1} A T$.
    B : sp.Matrix
        Canonical input matrix $B_{\text{can}} = T^{-1} B$.
    C : sp.Matrix
        Canonical output matrix $C_{\text{can}} = C T$.
    D : sp.Matrix
        Direct feedthrough matrix $D$.
    T : sp.Matrix | None
        Similarity transformation matrix $T$ mapping canonical states $z$ to original states $x = T z$.
    T_inv : sp.Matrix | None
        Inverse transformation matrix $T^{-1}$ mapping original states to canonical states $z = T^{-1} x$.
    form_type : str
        Type of canonical form (``'controllable'``, ``'observable'``, ``'jordan'``).
    is_valid : bool
        Whether the canonical transformation is nonsingular and valid.
    explanation : str
        Analytical derivation and pedagogical notes for this canonical form.
    dt : float | None, default=None
        Sampling time for discrete systems (None for continuous-time).
    """

    def __init__(
        self,
        A: sp.Matrix,
        B: sp.Matrix,
        C: sp.Matrix,
        D: sp.Matrix,
        T: sp.Matrix | None,
        T_inv: sp.Matrix | None,
        form_type: str,
        is_valid: bool = True,
        explanation: str = "",
        dt: float | None = None,
    ) -> None:
        self.A = A
        self.B = B
        self.C = C
        self.D = D
        self.T = T
        self.T_inv = T_inv
        self.form_type = form_type
        self.is_valid = is_valid
        self.explanation = explanation
        self.dt = dt

    def to_ss(self) -> StateSpace:
        """Convert the symbolic canonical form to a numerical StateSpace model.

        Returns
        -------
        StateSpace
            Numerical StateSpace instance.

        Raises
        ------
        ValueError
            If matrices contain unresolved symbolic variables.
        """
        from ctrlpy.models.state_space import StateSpace

        def _to_np(m: sp.Matrix) -> np.ndarray:
            arr = np.zeros((m.rows, m.cols), dtype=np.float64)
            for i in range(m.rows):
                for j in range(m.cols):
                    val = m[i, j]
                    if hasattr(val, "free_symbols") and val.free_symbols:
                        raise ValueError(
                            f"Cannot convert symbolic expression {val} to numeric float."
                        )
                    arr[i, j] = float(sp.re(val.evalf()))
            return arr

        return StateSpace(_to_np(self.A), _to_np(self.B), _to_np(self.C), _to_np(self.D))

    def _repr_latex_(self) -> str:
        r"""Render LaTeX representation for Jupyter environments."""
        name_map = {
            "controllable": "Controllable Canonical Form (Phase-Variable)",
            "observable": "Observable Canonical Form",
            "jordan": "Jordan / Diagonal Modal Canonical Form",
        }
        title = name_map.get(self.form_type, f"{self.form_type.title()} Canonical Form")
        a_tex = _sp_matrix_to_latex(self.A)
        b_tex = _sp_matrix_to_latex(self.B)
        c_tex = _sp_matrix_to_latex(self.C)
        d_tex = _sp_matrix_to_latex(self.D)

        t_tex = _sp_matrix_to_latex(self.T) if self.T is not None else r"\text{None (Singular)}"
        t_inv_tex = (
            _sp_matrix_to_latex(self.T_inv) if self.T_inv is not None else r"\text{None (Singular)}"
        )

        dt_sym = r"\dot{z}" if self.dt is None else "z[k+1]"
        state_sym = "z" if self.dt is None else "z[k]"
        u_sym = "u" if self.dt is None else "u[k]"
        y_sym = "y" if self.dt is None else "y[k]"

        lines = [
            r"\begin{aligned}",
            rf"\textbf{{{title}:}} \\",
            rf"{dt_sym} &= {a_tex} {state_sym} + {b_tex} {u_sym} \\",
            rf"{y_sym} &= {c_tex} {state_sym} + {d_tex} {u_sym} \\",
            rf"T &= {t_tex}, \quad T^{{-1}} = {t_inv_tex}",
            r"\end{aligned}",
        ]
        return "\n".join(lines)

    def _repr_markdown_(self) -> str:
        """Render Markdown for Jupyter environments."""
        return f"$${self._repr_latex_()}$$"

    def __str__(self) -> str:
        """ASCII representation of canonical form matrices."""
        return (
            f"=== {self.form_type.title()} Canonical Form ===\n"
            f"A =\n{self.A}\n\n"
            f"B =\n{self.B}\n\n"
            f"C =\n{self.C}\n\n"
            f"D =\n{self.D}\n\n"
            f"Transformation Matrix T (x = T z):\n{self.T}\n\n"
            f"Inverse Transformation T_inv (z = T^-1 x):\n{self.T_inv}"
        )

    def __repr__(self) -> str:
        return f"CanonicalFormResult(form_type={self.form_type!r}, is_valid={self.is_valid})"

to_ss

to_ss()

Convert the symbolic canonical form to a numerical StateSpace model.

Returns:

Type	Description
StateSpace	Numerical StateSpace instance.

Raises:

Type	Description
ValueError	If matrices contain unresolved symbolic variables.

Source code in src/ctrlpy/symbolic/state_space.py

def to_ss(self) -> StateSpace:
    """Convert the symbolic canonical form to a numerical StateSpace model.

    Returns
    -------
    StateSpace
        Numerical StateSpace instance.

    Raises
    ------
    ValueError
        If matrices contain unresolved symbolic variables.
    """
    from ctrlpy.models.state_space import StateSpace

    def _to_np(m: sp.Matrix) -> np.ndarray:
        arr = np.zeros((m.rows, m.cols), dtype=np.float64)
        for i in range(m.rows):
            for j in range(m.cols):
                val = m[i, j]
                if hasattr(val, "free_symbols") and val.free_symbols:
                    raise ValueError(
                        f"Cannot convert symbolic expression {val} to numeric float."
                    )
                arr[i, j] = float(sp.re(val.evalf()))
        return arr

    return StateSpace(_to_np(self.A), _to_np(self.B), _to_np(self.C), _to_np(self.D))

__str__

__str__()

ASCII representation of canonical form matrices.

Source code in src/ctrlpy/symbolic/state_space.py

def __str__(self) -> str:
    """ASCII representation of canonical form matrices."""
    return (
        f"=== {self.form_type.title()} Canonical Form ===\n"
        f"A =\n{self.A}\n\n"
        f"B =\n{self.B}\n\n"
        f"C =\n{self.C}\n\n"
        f"D =\n{self.D}\n\n"
        f"Transformation Matrix T (x = T z):\n{self.T}\n\n"
        f"Inverse Transformation T_inv (z = T^-1 x):\n{self.T_inv}"
    )

ModeAnalysis dataclass

Modal controllability and observability analysis for an individual eigenvalue.

Attributes:

Name	Type	Description
eigenvalue	Expr	The eigenvalue / pole \(\lambda_i\).
multiplicity	int	Algebraic multiplicity of \(\lambda_i\).
is_controllable	bool	True if \(\lambda_i\) satisfies the PBH controllability rank test.
is_observable	bool	True if \(\lambda_i\) satisfies the PBH observability rank test.
pbh_c_rank	int	Rank of the PBH controllability matrix \([\lambda_i I - A \quad B]\).
pbh_o_rank	int	Rank of the PBH observability matrix \(\begin{bmatrix} \lambda_i I - A \\ C \end{bmatrix}\).
kalman_type	str	Kalman decomposition classification: - 'co': Controllable and Observable - 'c_unobs': Controllable and Unobservable - 'unctrl_obs': Uncontrollable and Observable - 'unctrl_unobs': Uncontrollable and Unobservable
description	str	Human-readable explanation of the mode's physical behavior and pole-zero cancellation status.

Source code in src/ctrlpy/symbolic/state_space.py

@dataclass
class ModeAnalysis:
    r"""Modal controllability and observability analysis for an individual eigenvalue.

    Attributes
    ----------
    eigenvalue : sp.Expr
        The eigenvalue / pole $\lambda_i$.
    multiplicity : int
        Algebraic multiplicity of $\lambda_i$.
    is_controllable : bool
        True if $\lambda_i$ satisfies the PBH controllability rank test.
    is_observable : bool
        True if $\lambda_i$ satisfies the PBH observability rank test.
    pbh_c_rank : int
        Rank of the PBH controllability matrix $[\lambda_i I - A \quad B]$.
    pbh_o_rank : int
        Rank of the PBH observability matrix $\begin{bmatrix} \lambda_i I - A \\ C \end{bmatrix}$.
    kalman_type : str
        Kalman decomposition classification:
        - ``'co'``: Controllable and Observable
        - ``'c_unobs'``: Controllable and Unobservable
        - ``'unctrl_obs'``: Uncontrollable and Observable
        - ``'unctrl_unobs'``: Uncontrollable and Unobservable
    description : str
        Human-readable explanation of the mode's physical behavior and pole-zero cancellation status.
    """

    eigenvalue: sp.Expr
    multiplicity: int
    is_controllable: bool
    is_observable: bool
    pbh_c_rank: int
    pbh_o_rank: int
    kalman_type: str
    description: str

state_space_tutor

state_space_tutor(
    sys_or_A=None,
    B=None,
    C=None,
    D=None,
    *,
    num=None,
    den=None,
    dt=None,
)

Create a pedagogical StateSpaceTutor for canonical transformations and controllability/observability analysis.

Parameters:

Name	Type	Description	Default
sys_or_A	Any	StateSpace model, TransferFunction model, or \(A\) matrix.	None
B	Any	Input matrix \(B\).	None
C	Any	Output matrix \(C\).	None
D	Any	Feedthrough matrix \(D\).	None
num	Sequence[Any] | None	Numerator polynomial coefficients.	None
den	Sequence[Any] | None	Denominator polynomial coefficients.	None
dt	float | None	Discrete sampling period (None for continuous-time).	None

Returns:

Type	Description
StateSpaceTutor	The analyzed tutor instance.

Source code in src/ctrlpy/symbolic/state_space.py

def state_space_tutor(
    sys_or_A: Any = None,
    B: Any = None,
    C: Any = None,
    D: Any = None,
    *,
    num: Sequence[Any] | None = None,
    den: Sequence[Any] | None = None,
    dt: float | None = None,
) -> StateSpaceTutor:
    r"""Create a pedagogical StateSpaceTutor for canonical transformations and controllability/observability analysis.

    Parameters
    ----------
    sys_or_A : Any, optional
        StateSpace model, TransferFunction model, or $A$ matrix.
    B : Any, optional
        Input matrix $B$.
    C : Any, optional
        Output matrix $C$.
    D : Any, optional
        Feedthrough matrix $D$.
    num : Sequence[Any] | None, default=None
        Numerator polynomial coefficients.
    den : Sequence[Any] | None, default=None
        Denominator polynomial coefficients.
    dt : float | None, default=None
        Discrete sampling period (None for continuous-time).

    Returns
    -------
    StateSpaceTutor
        The analyzed tutor instance.
    """
    return StateSpaceTutor(sys_or_A, B, C, D, num=num, den=den, dt=dt)

controllable_canonical_form

controllable_canonical_form(
    sys_or_A=None,
    B=None,
    C=None,
    D=None,
    *,
    num=None,
    den=None,
    dt=None,
)

Compute the Controllable Canonical Form (Phase-Variable Form) of a system.

Parameters:

Name	Type	Description	Default
sys_or_A	Any	StateSpace model, TransferFunction model, or \(A\) matrix.	None
B	Any	Input matrix \(B\).	None
C	Any	Output matrix \(C\).	None
D	Any	Feedthrough matrix \(D\).	None
num	Sequence[Any] | None	Numerator polynomial coefficients.	None
den	Sequence[Any] | None	Denominator polynomial coefficients.	None
dt	float | None	Sampling period.	None

Returns:

Type	Description
CanonicalFormResult	Controllable canonical form representation with transformation matrix \(T_c\).

Source code in src/ctrlpy/symbolic/state_space.py

def controllable_canonical_form(
    sys_or_A: Any = None,
    B: Any = None,
    C: Any = None,
    D: Any = None,
    *,
    num: Sequence[Any] | None = None,
    den: Sequence[Any] | None = None,
    dt: float | None = None,
) -> CanonicalFormResult:
    """Compute the Controllable Canonical Form (Phase-Variable Form) of a system.

    Parameters
    ----------
    sys_or_A : Any, optional
        StateSpace model, TransferFunction model, or $A$ matrix.
    B : Any, optional
        Input matrix $B$.
    C : Any, optional
        Output matrix $C$.
    D : Any, optional
        Feedthrough matrix $D$.
    num : Sequence[Any] | None, default=None
        Numerator polynomial coefficients.
    den : Sequence[Any] | None, default=None
        Denominator polynomial coefficients.
    dt : float | None, default=None
        Sampling period.

    Returns
    -------
    CanonicalFormResult
        Controllable canonical form representation with transformation matrix $T_c$.
    """
    tutor = StateSpaceTutor(sys_or_A, B, C, D, num=num, den=den, dt=dt)
    return tutor.controllable_canonical_form()

observable_canonical_form

observable_canonical_form(
    sys_or_A=None,
    B=None,
    C=None,
    D=None,
    *,
    num=None,
    den=None,
    dt=None,
)

Compute the Observable Canonical Form of a system.

Parameters:

Name	Type	Description	Default
sys_or_A	Any	StateSpace model, TransferFunction model, or \(A\) matrix.	None
B	Any	Input matrix \(B\).	None
C	Any	Output matrix \(C\).	None
D	Any	Feedthrough matrix \(D\).	None
num	Sequence[Any] | None	Numerator polynomial coefficients.	None
den	Sequence[Any] | None	Denominator polynomial coefficients.	None
dt	float | None	Sampling period.	None

Returns:

Type	Description
CanonicalFormResult	Observable canonical form representation with transformation matrix \(T_o\).

Source code in src/ctrlpy/symbolic/state_space.py

def observable_canonical_form(
    sys_or_A: Any = None,
    B: Any = None,
    C: Any = None,
    D: Any = None,
    *,
    num: Sequence[Any] | None = None,
    den: Sequence[Any] | None = None,
    dt: float | None = None,
) -> CanonicalFormResult:
    """Compute the Observable Canonical Form of a system.

    Parameters
    ----------
    sys_or_A : Any, optional
        StateSpace model, TransferFunction model, or $A$ matrix.
    B : Any, optional
        Input matrix $B$.
    C : Any, optional
        Output matrix $C$.
    D : Any, optional
        Feedthrough matrix $D$.
    num : Sequence[Any] | None, default=None
        Numerator polynomial coefficients.
    den : Sequence[Any] | None, default=None
        Denominator polynomial coefficients.
    dt : float | None, default=None
        Sampling period.

    Returns
    -------
    CanonicalFormResult
        Observable canonical form representation with transformation matrix $T_o$.
    """
    tutor = StateSpaceTutor(sys_or_A, B, C, D, num=num, den=den, dt=dt)
    return tutor.observable_canonical_form()

jordan_canonical_form

jordan_canonical_form(
    sys_or_A=None,
    B=None,
    C=None,
    D=None,
    *,
    num=None,
    den=None,
    dt=None,
)

Compute the Jordan / Diagonal Modal Form of a system.

Parameters:

Name	Type	Description	Default
sys_or_A	Any	StateSpace model, TransferFunction model, or \(A\) matrix.	None
B	Any	Input matrix \(B\).	None
C	Any	Output matrix \(C\).	None
D	Any	Feedthrough matrix \(D\).	None
num	Sequence[Any] | None	Numerator polynomial coefficients.	None
den	Sequence[Any] | None	Denominator polynomial coefficients.	None
dt	float | None	Sampling period.	None

Returns:

Type	Description
CanonicalFormResult	Jordan / diagonal canonical form representation with modal matrix \(V\).

Source code in src/ctrlpy/symbolic/state_space.py

def jordan_canonical_form(
    sys_or_A: Any = None,
    B: Any = None,
    C: Any = None,
    D: Any = None,
    *,
    num: Sequence[Any] | None = None,
    den: Sequence[Any] | None = None,
    dt: float | None = None,
) -> CanonicalFormResult:
    """Compute the Jordan / Diagonal Modal Form of a system.

    Parameters
    ----------
    sys_or_A : Any, optional
        StateSpace model, TransferFunction model, or $A$ matrix.
    B : Any, optional
        Input matrix $B$.
    C : Any, optional
        Output matrix $C$.
    D : Any, optional
        Feedthrough matrix $D$.
    num : Sequence[Any] | None, default=None
        Numerator polynomial coefficients.
    den : Sequence[Any] | None, default=None
        Denominator polynomial coefficients.
    dt : float | None, default=None
        Sampling period.

    Returns
    -------
    CanonicalFormResult
        Jordan / diagonal canonical form representation with modal matrix $V$.
    """
    tutor = StateSpaceTutor(sys_or_A, B, C, D, num=num, den=den, dt=dt)
    return tutor.jordan_canonical_form()

controllability_matrix

controllability_matrix(sys_or_A=None, B=None)

Compute the symbolic Kalman controllability matrix \(\mathcal{C} = [B \quad AB \dots A^{n-1}B]\).

Parameters:

Name	Type	Description	Default
sys_or_A	Any	StateSpace model, or state transition matrix \(A\).	None
B	Any	Input matrix \(B\).	None

Returns:

Type	Description
Matrix	Controllability matrix \(\mathcal{C}\).

Source code in src/ctrlpy/symbolic/state_space.py

def controllability_matrix(
    sys_or_A: Any = None,
    B: Any = None,
) -> sp.Matrix:
    r"""Compute the symbolic Kalman controllability matrix $\mathcal{C} = [B \quad AB \dots A^{n-1}B]$.

    Parameters
    ----------
    sys_or_A : Any
        StateSpace model, or state transition matrix $A$.
    B : Any, optional
        Input matrix $B$.

    Returns
    -------
    sp.Matrix
        Controllability matrix $\mathcal{C}$.
    """
    if hasattr(sys_or_A, "A") and hasattr(sys_or_A, "B"):
        A_sp = _to_sympy_matrix(sys_or_A.A)
        B_sp = _to_sympy_matrix(sys_or_A.B)
    else:
        A_sp = _to_sympy_matrix(sys_or_A)
        if B is None:
            raise ValueError("Input matrix B must be provided when A is given.")
        B_sp = _to_sympy_matrix(B)

    n = A_sp.rows
    cols: list[sp.Matrix] = []
    for k in range(n):
        cols.append((A_sp**k) * B_sp)
    return sp.Matrix.hstack(*cols)

observability_matrix

observability_matrix(sys_or_A=None, C=None)

Compute the symbolic Kalman observability matrix \(\mathcal{O} = [C^T \quad A^TC^T \dots (A^{n-1})^TC^T]^T\).

Parameters:

Name	Type	Description	Default
sys_or_A	Any	StateSpace model, or state transition matrix \(A\).	None
C	Any	Output matrix \(C\).	None

Returns:

Type	Description
Matrix	Observability matrix \(\mathcal{O}\).

Source code in src/ctrlpy/symbolic/state_space.py

def observability_matrix(
    sys_or_A: Any = None,
    C: Any = None,
) -> sp.Matrix:
    r"""Compute the symbolic Kalman observability matrix $\mathcal{O} = [C^T \quad A^TC^T \dots (A^{n-1})^TC^T]^T$.

    Parameters
    ----------
    sys_or_A : Any
        StateSpace model, or state transition matrix $A$.
    C : Any, optional
        Output matrix $C$.

    Returns
    -------
    sp.Matrix
        Observability matrix $\mathcal{O}$.
    """
    if hasattr(sys_or_A, "A") and hasattr(sys_or_A, "C"):
        A_sp = _to_sympy_matrix(sys_or_A.A)
        C_sp = _to_sympy_matrix(sys_or_A.C)
    else:
        A_sp = _to_sympy_matrix(sys_or_A)
        if C is None:
            raise ValueError("Output matrix C must be provided when A is given.")
        C_sp = _to_sympy_matrix(C)

    n = A_sp.rows
    rows: list[sp.Matrix] = []
    for k in range(n):
        rows.append(C_sp * (A_sp**k))
    return sp.Matrix.vstack(*rows)