diff --git a/README.md b/README.md index 5abdd88..638adbe 100644 --- a/README.md +++ b/README.md @@ -27,6 +27,7 @@ PathSim-RF extends the [PathSim](https://github.com/pathsim/pathsim) simulation | Block | Description | Key Parameters | |-------|-------------|----------------| | `RFNetwork` | N-port network from S-parameter data (Touchstone) via vector fitting | `ntwk`, `auto_fit` | +| `NetlistStateSpace` | Conversion of spice netlist into pathsim StateSpace | `netlist`, `output_voltages`, `output_currents` | | `TransmissionLine` | Lossy delay-based transmission line (scattering domain) | `length`, `er`, `attenuation`, `Z0` | | `RFAmplifier` | Amplifier with optional IP3 nonlinearity | `gain` [dB], `IIP3` [dBm], `P1dB` [dBm], `Z0` | | `RFMixer` | Ideal frequency converter (time-domain multiplication) | `conversion_gain` [dB], `Z0` | diff --git a/docs/source/examples/netlist_import.ipynb b/docs/source/examples/netlist_import.ipynb new file mode 100644 index 0000000..35fee01 --- /dev/null +++ b/docs/source/examples/netlist_import.ipynb @@ -0,0 +1,301 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Netlist State-Space Validation\n", + "\n", + "This example validates the `NetlistStateSpace` block against an equivalent PathSim implementation of an LC filter with a resistive load.\n", + "\n", + "Two circuit descriptions are used. The first netlist contains the load explicitly as a resistor. The second represents the same load using a behavioral current source. Both netlists are passed directly to `NetlistStateSpace` without requiring external `.net` files.\n", + "\n", + "The simulation compares the PathSim implementation and the imported netlist models using the recorded capacitor voltage and load current signals.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "from pathsim import Simulation, Connection\n", + "from pathsim.blocks import (\n", + " Adder,\n", + " Amplifier,\n", + " Integrator,\n", + " PulseSource,\n", + " Scope\n", + ")\n", + "from pathsim.solvers import EUB\n", + "from pathsim_rf.netlist_to_statespace import NetlistStateSpace\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Rload = 50\n", + "Cfilter = 50e-6\n", + "Lfilter = 0.5e-3\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Netlists\n", + "\n", + "The first netlist contains the 50 Ω load explicitly as `Rload`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "filter_netlist = \"\"\"\n", + "*validation netlist for pathsim\n", + ".param FREQ = 9k\n", + ".param DUTY = 0.5\n", + ".param V_LOW = 0\n", + ".param V_HIGH = 1\n", + ".param T_PER = {1/FREQ}\n", + ".param T_ON = {T_PER * DUTY}\n", + ".param T_RISE = 10n\n", + ".param T_FALL = 10n\n", + "\n", + "V1 in 0 PULSE({V_LOW} {V_HIGH} 0 {T_RISE} {T_FALL} {T_ON} {T_PER})\n", + "Lfilter in n1 0.5e-3\n", + "Cfilter n1 0 50e-6\n", + "Rload n1 0 50\n", + ".tran 0.00001 0.2 0 0.000001\n", + ".end\n", + "\"\"\"\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The second netlist represents the load using a behavioral current source.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "filter_without_load_netlist = \"\"\"\n", + "*validation netlist for pathsim\n", + ".param FREQ = 9k\n", + ".param DUTY = 0.5\n", + ".param V_LOW = 0\n", + ".param V_HIGH = 1\n", + ".param T_PER = {1/FREQ}\n", + ".param T_ON = {T_PER * DUTY}\n", + ".param T_RISE = 10n\n", + ".param T_FALL = 10n\n", + "\n", + "V1 in 0 PULSE({V_LOW} {V_HIGH} 0 {T_RISE} {T_FALL} {T_ON} {T_PER})\n", + "Bload n1 0 I=V(n1)/50\n", + "Lfilter in n1 0.5e-3\n", + "Cfilter n1 0 50e-6\n", + ".tran 0.00001 0.2 0 0.00001\n", + ".end\n", + "\"\"\"\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Blocks\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Sources\n", + "pulsesource = PulseSource(\n", + " T=1/9000,\n", + " amplitude=1,\n", + " duty=0.5,\n", + " t_rise=10e-9, \n", + " t_fall=10e-9\n", + ")\n", + "\n", + "state_space = NetlistStateSpace(\n", + " filter_netlist,\n", + " output_voltages=[\"n1\"],\n", + " output_currents=[\"Rload\"],\n", + ")\n", + "\n", + "state_space_bis = NetlistStateSpace(\n", + " filter_without_load_netlist,\n", + " output_voltages=[\"n1\"],\n", + " output_currents=[\"Lfilter\"],\n", + ")\n", + "\n", + "# Dynamic\n", + "integrator = Integrator(initial_value=0)\n", + "block_2 = Integrator()\n", + "\n", + "# Algebraic\n", + "adder = Adder(operations=\"+-\")\n", + "amplifier = Amplifier(gain=1/Lfilter)\n", + "block_5 = Adder(operations=\"+-\")\n", + "block_6 = Amplifier(gain=1/Cfilter)\n", + "block_7 = Amplifier(gain=1/Rload)\n", + "block_8 = Amplifier(gain=1/Rload)\n", + "\n", + "# Recording\n", + "scope = Scope(\n", + " sampling_period=1/100000,\n", + " labels=[\n", + " 'V(pwm)',\n", + " 'V(Cfull_pathsim)',\n", + " 'V(C_net_imported)',\n", + " 'V(Bload)'\n", + " ]\n", + ")\n", + "\n", + "blocks = [\n", + " pulsesource,\n", + " integrator,\n", + " block_2,\n", + " adder,\n", + " amplifier,\n", + " block_5,\n", + " block_6,\n", + " block_7,\n", + " scope,\n", + " state_space,\n", + " state_space_bis,\n", + " block_8\n", + "]\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Connections\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "conn_0 = Connection(pulsesource[0], adder[0], state_space[0], state_space_bis[0])\n", + "conn_1 = Connection(adder[0], integrator[0])\n", + "conn_2 = Connection(integrator[0], amplifier[0])\n", + "conn_3 = Connection(amplifier[0], block_5[0])\n", + "conn_4 = Connection(block_5[0], block_2[0])\n", + "conn_5 = Connection(block_2[0], block_6[0])\n", + "conn_6 = Connection(block_6[0], block_7[0])\n", + "conn_7 = Connection(block_7[0], block_5[1])\n", + "conn_8 = Connection(block_6[0], adder[1])\n", + "conn_9 = Connection(pulsesource[0], scope[0])\n", + "conn_10 = Connection(block_6[0], scope[1])\n", + "conn_11 = Connection(state_space[0], scope[2])\n", + "conn_12 = Connection(state_space_bis[0], block_8[0], scope[3])\n", + "conn_13 = Connection(block_8[0], state_space_bis[1])\n", + "\n", + "connections = [\n", + " conn_0,\n", + " conn_1,\n", + " conn_2,\n", + " conn_3,\n", + " conn_4,\n", + " conn_5,\n", + " conn_6,\n", + " conn_7,\n", + " conn_8,\n", + " conn_9,\n", + " conn_10,\n", + " conn_11,\n", + " conn_12,\n", + " conn_13\n", + "]\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Simulation\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sim = Simulation(\n", + " blocks,\n", + " connections,\n", + " Solver=EUB,\n", + " dt=0.5e-6,\n", + " dt_min=1e-16,\n", + " dt_max=1e-6,\n", + " tolerance_lte_rel=1e-6,\n", + " tolerance_lte_abs=1e-10,\n", + " tolerance_fpi=1e-10,\n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Results\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sim.run(duration=0.05)\n", + "\n", + "sim.plot()\n", + "plt.show()\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.x", + "mimetype": "text/x-python", + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "pygments_lexer": "ipython3", + "nbconvert_exporter": "python", + "file_extension": ".py" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/pathsim_rf/__init__.py b/src/pathsim_rf/__init__.py index 0531686..c3fe14b 100644 --- a/src/pathsim_rf/__init__.py +++ b/src/pathsim_rf/__init__.py @@ -15,6 +15,7 @@ from .transmission_line import * from .amplifier import * from .mixer import * +from .netlist_to_statespace import NetlistStateSpace try: from .network import * diff --git a/src/pathsim_rf/netlist_to_statespace.py b/src/pathsim_rf/netlist_to_statespace.py new file mode 100644 index 0000000..0cb3743 --- /dev/null +++ b/src/pathsim_rf/netlist_to_statespace.py @@ -0,0 +1,831 @@ +""" +netlist_to_statespace.py + +Turns a small SPICE-like netlist of R, L, C, K (mutual inductance / coupling +coefficient), V (independent voltage source) and I (independent current +source) elements into an explicit LTI state-space model + + xdot = A x + B u + y = C x + D u + +via Modified Nodal Analysis (MNA) + elimination of the algebraic +(non-storage) unknowns. This is the standard "state-variable method" used +to derive circuit ODEs by hand, automated in code so it scales to +messy topologies and handles coupled inductors correctly. + +Only LINEAR elements are supported (no diodes/switches/etc - PathSim itself +handles those fine via its algebraic-loop solver and Function blocks, see +the accompanying explanation). +""" + +from __future__ import annotations +import numpy as np +from dataclasses import dataclass +from pathlib import Path +import re +import warnings + +from pathsim.blocks.lti import StateSpace + + +# -------------------------------------------------------------------------- +# Netlist parsing +# -------------------------------------------------------------------------- + +@dataclass +class Element: + kind: str # 'R','L','C','V','I','K' + name: str + n1: str = None + n2: str = None + value: float = None + # for K elements: + l1: str = None + l2: str = None + + +SI_MULTIPLIERS = { + "f": 1e-15, + "p": 1e-12, + "n": 1e-9, + "u": 1e-6, + "µ": 1e-6, + "m": 1e-3, + "": 1.0, + "k": 1e3, + "K": 1e3, + "meg": 1e6, + "Meg": 1e6, + "M": 1e6, + "g": 1e9, + "G": 1e9, + "t": 1e12, + "T": 1e12, +} + +_PARAM_PATTERN = re.compile(r"^\.param\s+([A-Za-z_]\w*)\s*=\s*(.+)$") +_VALUE_PATTERN = re.compile( + r"^\s*([+-]?\d*\.?\d+(?:[eE][+-]?\d+)?)([A-Za-zµ]*)\s*$" +) +_VALID_KINDS = {"R", "L", "C", "V", "I", "K"} +_IGNORED_LINE_STARTS = (".", "*", '"') + + +def _normalize_node(node: str) -> str: + """Normalize ground aliases to canonical node '0'.""" + return "0" if node.strip().lower() in {"0", "gnd"} else node + + +def _get_si_multiplier(suffix: str) -> float: + """Return numeric multiplier for an SI prefix/suffix string.""" + if suffix in SI_MULTIPLIERS: + return SI_MULTIPLIERS[suffix] + if suffix.lower() in SI_MULTIPLIERS: + return SI_MULTIPLIERS[suffix.lower()] + + candidates = sorted(SI_MULTIPLIERS.keys(), key=len, reverse=True) + for key in candidates: + if not key: + continue + if suffix.startswith(key): + return SI_MULTIPLIERS[key] + low_key = key.lower() + if suffix.lower().startswith(low_key): + return SI_MULTIPLIERS[low_key] + raise ValueError(f"Unknown SI prefix: '{suffix}'") + + +def parse_value_with_units(value_str: str, params: dict[str, str] | None = None) -> float: + """ + Parse SPICE-like numeric tokens with SI prefixes and optional unit tails. + + Supported forms include: + - plain/scientific floats: ``1``, ``100.e-3``, ``2.2e6`` + - SI-prefixed tokens: ``4u``, ``10k``, ``3Meg``, ``2.2kOhm``, ``1mH`` + - parameter references: ``{RMAIN}`` when ``params`` is provided + """ + if value_str is None: + raise ValueError("Missing value") + + value = value_str.strip() + try: + return float(value) + except ValueError: + pass + + if value.startswith("{") and value.endswith("}"): + if params is None: + raise ValueError(f"Unknown parameter reference: '{value}'") + name = value[1:-1].strip() + if name not in params: + raise ValueError(f"Unknown parameter: '{name}'") + return parse_value_with_units(params[name], params) + + match = _VALUE_PATTERN.fullmatch(value) + if match is None: + raise ValueError(f"Invalid value format: '{value_str}'") + + number, suffix = match.groups() + return float(number) * _get_si_multiplier(suffix) + + +def _parse_optional_source_value(raw: str, params: dict[str, str]) -> float | None: + """Parse source literal values, returning None for waveform expressions.""" + try: + return parse_value_with_units(raw, params) + except ValueError: + # Many source statements use waveforms (PWL, SIN, EXP, etc.); those + # are runtime waveforms and do not affect linearization here. + return None + + +def _parse_behavioral_source(parts: list[str], line_number: int, line: str) -> Element: + """ + Parse LTspice behavioral source shorthand: + Bx n+ n- I= -> modeled as current source input placeholder + Bx n+ n- V= -> modeled as voltage source input placeholder + """ + if len(parts) < 4: + raise ValueError(f"Malformed behavioral source at line {line_number}: '{line}'") + + name = parts[0] + n1, n2 = _normalize_node(parts[1]), _normalize_node(parts[2]) + expr = " ".join(parts[3:]).strip() + expr_upper = expr.upper() + if expr_upper.startswith("I="): + return Element(kind="I", name=name, n1=n1, n2=n2, value=None) + if expr_upper.startswith("V="): + return Element(kind="V", name=name, n1=n1, n2=n2, value=None) + raise ValueError( + f"Unsupported behavioral source expression at line {line_number}: '{line}'" + ) + + +def parse_netlist(text: str) -> list[Element]: + """ + Parse netlist text into linear element records. + + One element per line, whitespace separated; '#' or ';' starts a comment. + R n1 n2 value + L n1 n2 value + C n1 n2 value + V n+ n- value (value is a placeholder; real waveform comes + from the PathSim Source block at sim time) + I n+ n- value (same) + B n+ n- I= (behavioral current source placeholder) + B n+ n- V= (behavioral voltage source placeholder) + K Lname1 Lname2 k (coupling coefficient, -1<=k<=1) + Node '0' (or 'gnd') is ground. + Deck/meta lines beginning with '.', '*', or a quoted header are ignored. + """ + params: dict[str, str] = {} + raw_element_lines: list[tuple[int, str]] = [] + + for line_number, raw in enumerate(text.splitlines(), start=1): + line = raw.split("#")[0].split(";")[0].strip() + if not line: + continue + match = _PARAM_PATTERN.fullmatch(line) + if match is not None: + params[match.group(1)] = match.group(2).strip() + continue + raw_element_lines.append((line_number, line)) + + elements: list[Element] = [] + for line_number, line in raw_element_lines: + parts = line.split() + if not parts: + continue + if parts[0].startswith(_IGNORED_LINE_STARTS): + continue + + name = parts[0] + kind = name[0].upper() + if kind == "B": + elements.append(_parse_behavioral_source(parts, line_number, line)) + continue + + if kind not in _VALID_KINDS: + raise ValueError( + f"Unsupported element '{name}' at line {line_number}: '{line}'" + ) + + if kind == 'K': + if len(parts) < 4: + raise ValueError(f"Malformed K element at line {line_number}: '{line}'") + elements.append(Element(kind='K', name=name, l1=parts[1], l2=parts[2], + value=parse_value_with_units(parts[3], params))) + else: + if len(parts) < 3: + raise ValueError(f"Malformed element at line {line_number}: '{line}'") + n1, n2 = _normalize_node(parts[1]), _normalize_node(parts[2]) + val_token = " ".join(parts[3:]) if len(parts) > 3 else None + if kind in {"R", "L", "C"}: + if val_token is None: + raise ValueError(f"Missing value for element '{name}' at line {line_number}") + val = parse_value_with_units(val_token, params) + else: + val = _parse_optional_source_value(val_token, params) if val_token else None + elements.append(Element(kind=kind, name=name, n1=n1, n2=n2, value=val)) + return elements + + +def parse_netlist_file(path: str | Path) -> list[Element]: + """Load and parse a netlist file.""" + with open(path, "r", encoding="utf-8") as handle: + return parse_netlist(handle.read()) + + +GND = {'0', 'gnd', 'GND'} + + +# -------------------------------------------------------------------------- +# MNA DAE -> explicit state space reduction +# -------------------------------------------------------------------------- + +class CircuitModel: + """ + Build an explicit state-space model from a linear circuit netlist. + + The circuit is first stamped as the modified nodal analysis (MNA) system + + ``E * dz/dt + G * z = B * u``, + + where ``z`` contains node voltages, inductor currents, and ideal voltage + source currents. The MNA unknowns are partitioned into differential states + ``x_d`` and algebraic unknowns ``x_a``. Eliminating ``x_a`` gives + + ``x_a = Phi * x_d + Psi * u`` + + and the final explicit model + + ``dx_d/dt = A * x_d + B_ss * u``. + + The MNA matrices are stamped directly into floating-point arrays and + reduced using SciPy sparse LU solves. Singular solves are retried with + bounded diagonal regularization and emit :class:`RuntimeWarning` when + regularization is used. + + Parameters + ---------- + elements: + Output of :func:`parse_netlist` / :func:`parse_netlist_file`. + + Attributes + ---------- + A: + State matrix of the reduced explicit system. + B_ss: + Input matrix of the reduced explicit system. + Phi: + Map from differential states to eliminated algebraic variables. + Psi: + Map from inputs to eliminated algebraic variables. + state_labels: + Labels matching the row and column ordering of the state matrices. + """ + + def __init__(self, elements: list[Element]): + + self.elements = elements + + self.R = [e for e in elements if e.kind == 'R'] + self.L = [e for e in elements if e.kind == 'L'] + self.C = [e for e in elements if e.kind == 'C'] + self.V = [e for e in elements if e.kind == 'V'] + self.I = [e for e in elements if e.kind == 'I'] + self.K = [e for e in elements if e.kind == 'K'] + self.dipoles_by_name = { + e.name: e for e in elements if e.kind in {"R", "L", "C", "V", "I"} + } + + # ---- nodes ---- + nodes = [] + for e in elements: + if e.kind == 'K': + continue + for n in (e.n1, e.n2): + if n not in GND and n not in nodes: + nodes.append(n) + self.nodes = nodes # ordered list of non-ground node names + self.node_idx = {n: i for i, n in enumerate(nodes)} + self.n_nodes = len(nodes) + + self.L_names = [e.name for e in self.L] + self.L_idx = {name: i for i, name in enumerate(self.L_names)} + self.n_L = len(self.L) + + self.V_names = [e.name for e in self.V] + self.V_idx = {name: i for i, name in enumerate(self.V_names)} + self.n_V = len(self.V) + + # unknown ordering: [node voltages] + [inductor currents] + [V-source currents] + self.n_total = self.n_nodes + self.n_L + self.n_V + + def col_node(n): + return None if n in GND else self.node_idx[n] + + def col_iL(name): + return self.n_nodes + self.L_idx[name] + + def col_iV(name): + return self.n_nodes + self.n_L + self.V_idx[name] + + self._col_node, self._col_iL, self._col_iV = col_node, col_iL, col_iV + + # ---- inductance matrix (with mutual terms) ---- + Lmat = np.zeros((self.n_L, self.n_L), dtype=float) + l_values = {e.name: e.value for e in self.L} + for e in self.L: + i = self.L_idx[e.name] + Lmat[i, i] = float(e.value) + for k in self.K: + i, j = self.L_idx[k.l1], self.L_idx[k.l2] + Li = l_values[k.l1] + Lj = l_values[k.l2] + M = float(k.value) * np.sqrt(float(Li) * float(Lj)) + Lmat[i, j] += M + Lmat[j, i] += M + self.Lmat = Lmat + + # ---- inputs: one column per V source, then one column per I source ---- + self.input_names = self.V_names + [e.name for e in self.I] + self.n_u = len(self.input_names) + if self.n_u == 0: + raise ValueError( + "Circuit must contain at least one independent voltage or current source." + ) + + # ---- build E (dynamic) and G (algebraic) matrices, and B (input map) ---- + n = self.n_total + E = np.zeros((n, n), dtype=float) + G = np.zeros((n, n), dtype=float) + B = np.zeros((n, self.n_u), dtype=float) + + def stamp_G(row, col, val): + if row is not None and col is not None: + G[row, col] += val + + # Resistors: contribute to node KCL rows only + for e in self.R: + a, b = col_node(e.n1), col_node(e.n2) + g = 1.0 / float(e.value) + stamp_G(a, a, g); stamp_G(b, b, g) + stamp_G(a, b, -g); stamp_G(b, a, -g) + + # Capacitors: contribute dv/dt terms to node KCL rows (the E matrix) + for e in self.C: + a, b = col_node(e.n1), col_node(e.n2) + c = float(e.value) + if a is not None: E[a, a] += c + if b is not None: E[b, b] += c + if a is not None and b is not None: + E[a, b] -= c + E[b, a] -= c + + # Inductors: KCL stamp (current unknown enters/leaves nodes) + + # dedicated branch row v_na - v_nb - sum_j Lij * d(iLj)/dt = 0 + for e in self.L: + a, b = col_node(e.n1), col_node(e.n2) + iL_col = col_iL(e.name) + row = iL_col # branch row shares index with its current unknown + stamp_G(a, iL_col, 1) + stamp_G(b, iL_col, -1) + if a is not None: G[row, a] += 1 + if b is not None: G[row, b] -= 1 + # branch eqn: v_na - v_nb - L*d(iL)/dt - sum_j M_ij*d(iLj)/dt = 0 + # => E[row, iLj] = -Lij (note the minus sign!) + i = self.L_idx[e.name] + for j, name_j in enumerate(self.L_names): + Lij = self.Lmat[i, j] + if Lij != 0: + E[row, col_iL(name_j)] -= Lij + + # Voltage sources: KCL stamp + branch row v_na - v_nb = u(t) + for e in self.V: + a, b = col_node(e.n1), col_node(e.n2) + iV_col = col_iV(e.name) + row = iV_col + stamp_G(a, iV_col, 1) + stamp_G(b, iV_col, -1) + if a is not None: G[row, a] += 1 + if b is not None: G[row, b] -= 1 + u_col = self.input_names.index(e.name) + B[row, u_col] = 1 + + # Current sources: pure RHS injection into node KCL rows. + # Convention: positive I flows from n1 -> n2 *through the source*, + # i.e. it delivers current INTO n2 and draws it OUT of n1 from the + # external circuit's point of view. + for e in self.I: + a, b = col_node(e.n1), col_node(e.n2) + u_col = self.input_names.index(e.name) + if a is not None: B[a, u_col] -= 1 + if b is not None: B[b, u_col] += 1 + + self.E, self.G, self.B = E, G, B + + # ---- differential / algebraic partition ---- + diff_rows = [r for r in range(n) if any(E[r, c] != 0 for c in range(n))] + alg_rows = [r for r in range(n) if r not in diff_rows] + self.diff_rows, self.alg_rows = diff_rows, alg_rows + + # sanity: the state columns should be exactly node-voltages that own a + # nonzero E column plus all inductor currents; algebraic columns = rest + diff_cols = sorted(set(c for r in diff_rows for c in range(n) if E[r, c] != 0)) + alg_cols = [c for c in range(n) if c not in diff_cols] + if len(diff_cols) != len(diff_rows) or len(alg_cols) != len(alg_rows): + raise ValueError( + "Circuit is degenerate for this reduction (e.g. an all-capacitor " + "loop, an all-inductor cutset, or a floating node). " + f"diff_rows={len(diff_rows)} diff_cols={len(diff_cols)} " + f"alg_rows={len(alg_rows)} alg_cols={len(alg_cols)}" + ) + self.diff_cols, self.alg_cols = diff_cols, alg_cols + + self._set_state_labels() + self._reduce() + self._selected_output_rows: list[tuple[np.ndarray, np.ndarray]] = [] + self.output_labels: list[str] = [] + + def _set_state_labels(self): + """Set labels for the differential state vector in `self.diff_cols` order.""" + dc = self.diff_cols + labels = [] + for c in dc: + if c < self.n_nodes: + labels.append(f"v_{self.nodes[c]}") + else: + labels.append(f"i_{self.L_names[c - self.n_nodes]}") + self.state_labels = labels + self.state_label_idx = {label: i for i, label in enumerate(labels)} + + def _raise_algebraic_singular(self): + raise ValueError( + "Algebraic subsystem is singular. Usually means a loop made " + "purely of ideal voltage sources (and/or 0-ohm shorts), or a " + "node with no DC path to ground." + ) + + def _raise_state_singular(self): + raise ValueError( + "State/mass matrix is singular: this circuit's capacitor " + "voltages (or inductor currents) aren't independent, so the " + "naive 'one state per capacitor-touched node' selection " + "over-counted states. Classic cause: a capacitor whose *both* " + "terminals only reach the rest of the circuit through that " + "same capacitor (no other cap ties either node down " + "independently) -- e.g. 'R1 in a / C1 a b / R2 b 0' with " + "nothing else at a or b. Only one of v_a, v_b is really an " + "independent state there; the other is pinned by KCL. This " + "reduction doesn't do full tree/cotree state selection, so it " + "can't detect that automatically yet. Workarounds: (1) add a " + "TINY STRAY CAPACITANCE FROM ONE OF THE FLOATING NODES TO " + "GROUND (not a resistor -- the degeneracy lives in the " + "capacitor/mass matrix, a parallel resistor doesn't touch it " + "at all). e.g. 'Cstray b 0 1e-15' -- this gives that node's " + "own row independent rank; it adds one extra, extremely fast " + "eigenvalue (a numerical artifact, orders of magnitude faster " + "than your real dynamics) alongside the correct physical " + "pole(s), or (2) hand-pick the true independent capacitor " + "voltage as the state and eliminate the redundant node " + "yourself before building the netlist." + ) + + @staticmethod + def _safe_matmul(a: np.ndarray, b: np.ndarray) -> np.ndarray: + """ + Matrix multiply without relying on NumPy BLAS-backed matmul. + This avoids environment-specific native linear algebra crashes. + """ + if a.shape[1] != b.shape[0]: + raise ValueError(f"Incompatible shapes for matmul: {a.shape} and {b.shape}") + out = np.zeros((a.shape[0], b.shape[1]), dtype=float) + for i in range(a.shape[0]): + for k in range(a.shape[1]): + aik = a[i, k] + if aik != 0.0: + out[i, :] += aik * b[k, :] + return out + + @staticmethod + def _solve_with_scipy(mat: np.ndarray, rhs: np.ndarray, context: str) -> np.ndarray: + """Solve right-hand sides with SciPy sparse LU and bounded regularization.""" + import scipy.sparse as sps + from scipy.sparse.linalg import splu + + rhs_2d = rhs if rhs.ndim == 2 else rhs.reshape((-1, 1)) + if mat.shape[0] != mat.shape[1]: + raise ValueError(f"{context} matrix must be square, got {mat.shape}") + if rhs_2d.shape[0] != mat.shape[0]: + raise ValueError( + f"{context} rhs has incompatible shape {rhs_2d.shape} for matrix {mat.shape}" + ) + + last_error: Exception | None = None + for eps in (0.0, 1e-15, 1e-12, 1e-9, 1e-6): + try: + matrix = mat if eps == 0.0 else mat + eps * np.eye(mat.shape[0]) + factor = splu(sps.csc_matrix(matrix)) + columns = [ + np.asarray(factor.solve(np.asarray(column, dtype=float)), dtype=float) + for column in rhs_2d.T + ] + solution = np.column_stack(columns) + if eps > 0.0: + warnings.warn( + f"Circuit reduction regularized singular {context} with eps={eps:.1e}", + RuntimeWarning, + ) + return solution if rhs.ndim == 2 else solution[:, 0] + except (RuntimeError, ValueError) as exc: + last_error = exc + + raise ValueError(f"Failed to solve {context} with SciPy sparse LU.") from last_error + + def _reduce(self): + """Eliminate algebraic variables and solve the mass system with SciPy.""" + E, G, B = self.E, self.G, self.B + dr, ar = self.diff_rows, self.alg_rows + dc, ac = self.diff_cols, self.alg_cols + + E_dd = E[np.ix_(dr, dc)] + G_alg_d = G[np.ix_(ar, dc)] + G_alg_a = G[np.ix_(ar, ac)] + G_diff_d = G[np.ix_(dr, dc)] + G_diff_a = G[np.ix_(dr, ac)] + B_alg = B[ar, :] + B_diff = B[dr, :] + + try: + Phi = -self._solve_with_scipy(G_alg_a, G_alg_d, "algebraic subsystem") + Psi = self._solve_with_scipy(G_alg_a, B_alg, "algebraic subsystem") + except ValueError: + self._raise_algebraic_singular() + + try: + A = -self._solve_with_scipy( + E_dd, + G_diff_d + self._safe_matmul(G_diff_a, Phi), + "state/mass subsystem", + ) + Bmat = self._solve_with_scipy( + E_dd, + B_diff - self._safe_matmul(G_diff_a, Psi), + "state/mass subsystem", + ) + except ValueError: + self._raise_state_singular() + + self.A = A + self.B_ss = Bmat + self.Phi = Phi + self.Psi = Psi + + def _build_output_rows( + self, + coefficient_rows: list[dict[str, float]], + ) -> tuple[np.ndarray, np.ndarray]: + """ + Project raw MNA output expressions to state-space rows. + + Each coefficient mapping uses keys of the form ``node:name``, + ``iL:name``, or ``iV:name``. Algebraic unknowns are substituted through + ``Phi`` and ``Psi`` so each result satisfies ``y = C*x + D*u``. + """ + if not coefficient_rows: + return ( + np.empty((0, len(self.state_labels)), dtype=float), + np.empty((0, self.n_u), dtype=float), + ) + + vec = np.zeros((len(coefficient_rows), self.n_total), dtype=float) + for row, coeffs in enumerate(coefficient_rows): + for key, val in coeffs.items(): + try: + typ, name = key.split(":", 1) + except ValueError as exc: + raise ValueError(f"Invalid output coefficient key '{key}'") from exc + + if typ == "node": + col = self._col_node(name) + elif typ == "iL": + col = self._col_iL(name) + elif typ == "iV": + col = self._col_iV(name) + else: + raise ValueError(f"Unknown output coefficient type '{typ}'") + if col is not None: + vec[row, col] += float(val) + + vec_d = vec[:, self.diff_cols] + vec_a = vec[:, self.alg_cols] + C = vec_d + self._safe_matmul(vec_a, self.Phi) + D = self._safe_matmul(vec_a, self.Psi) + return C, D + + def _build_output_row( + self, + coefficients: dict[str, float], + ) -> tuple[np.ndarray, np.ndarray]: + """Project one raw MNA expression to one state-space output row.""" + C, D = self._build_output_rows([coefficients]) + return C[0], D[0] + + def _state_derivative_row( + self, + node: str, + ) -> tuple[np.ndarray, np.ndarray]: + """Return rows representing ``dV(node)/dt = C*x + D*u``.""" + n_states = len(self.state_labels) + if node in GND: + return np.zeros(n_states), np.zeros(self.n_u) + label = f"v_{node}" + if label not in self.state_labels: + raise ValueError( + f"node '{node}' has no capacitor attached, so its voltage isn't " + f"a state and dv/dt isn't defined this way. States: {self.state_labels}" + ) + idx = self.state_label_idx[label] + return self.A[idx, :], self.B_ss[idx, :] + + def add_node_voltage_output(self, node_name: str) -> None: + """ + Add a node-to-ground voltage to the selected system outputs. + + Parameters + ---------- + node_name: + Netlist node whose voltage is measured relative to ground. + + Raises + ------ + ValueError + If ``node_name`` is not present in the circuit. + """ + node = _normalize_node(node_name) + if node not in GND and node not in self.node_idx: + raise ValueError(f"Unknown node '{node_name}'") + + coefficients = {} if node in GND else {f"node:{node}": 1.0} + self._selected_output_rows.append(self._build_output_row(coefficients)) + self.output_labels.append(f"V({node})") + + def add_dipole_current_output(self, dipole_name: str) -> None: + """ + Add the current through a two-terminal netlist element. + + Parameters + ---------- + dipole_name: + Name of a resistor, inductor, capacitor, voltage source, or current + source. Positive current follows the element declaration from + ``n1`` to ``n2``. + + Raises + ------ + ValueError + If no supported dipole has the requested name. + + Notes + ----- + Resistor current is derived from Ohm's law. Inductor and voltage-source + currents are MNA branch unknowns. Capacitor current is computed from its + voltage derivative. Current-source current is its corresponding input + signal directly. + """ + try: + element = self.dipoles_by_name[dipole_name] + except KeyError as exc: + raise ValueError(f"Unknown dipole '{dipole_name}'") from exc + + if element.kind == "R": + coefficients = {} + if element.n1 not in GND: + coefficients[f"node:{element.n1}"] = 1.0 / element.value + if element.n2 not in GND: + key = f"node:{element.n2}" + coefficients[key] = coefficients.get(key, 0.0) - 1.0 / element.value + row = self._build_output_row(coefficients) + elif element.kind == "L": + row = self._build_output_row({f"iL:{element.name}": 1.0}) + elif element.kind == "C": + A_n1, B_n1 = self._state_derivative_row(element.n1) + A_n2, B_n2 = self._state_derivative_row(element.n2) + row = ( + element.value * (A_n1 - A_n2), + element.value * (B_n1 - B_n2), + ) + elif element.kind == "V": + row = self._build_output_row({f"iV:{element.name}": 1.0}) + elif element.kind == "I": + C_row = np.zeros(len(self.state_labels), dtype=float) + D_row = np.zeros(self.n_u, dtype=float) + D_row[self.input_names.index(element.name)] = 1.0 + row = C_row, D_row + else: + raise ValueError( + f"Element '{dipole_name}' of type '{element.kind}' is not a dipole" + ) + + self._selected_output_rows.append(row) + self.output_labels.append(f"I({element.name})") + + def get_system( + self, + ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """ + Return the reduced system with all selected outputs. + + Returns + ------- + A, B, C, D: + Matrices satisfying ``dx/dt = A*x + B*u`` and + ``y = C*x + D*u``. Rows of ``C`` and ``D`` follow the order in + which outputs were added. With no selected outputs, ``C`` and ``D`` + have zero rows and retain the correct state/input column counts. + """ + if not self._selected_output_rows: + return ( + self.A, + self.B_ss, + np.empty((0, len(self.state_labels)), dtype=float), + np.empty((0, self.n_u), dtype=float), + ) + C = np.vstack([row[0] for row in self._selected_output_rows]) + D = np.vstack([row[1] for row in self._selected_output_rows]) + return self.A, self.B_ss, C, D + + +class NetlistStateSpace(StateSpace): + """ + PathSim state-space block constructed directly from a linear netlist. + + Parameters + ---------- + netlist: + Existing netlist path or inline netlist text. A :class:`str` is treated + as a path when it names an existing file and as netlist text otherwise. + A :class:`Path` is always treated as an explicit file path. + output_voltages: + Node names exposed as node-to-ground voltage outputs. + output_currents: + Names of R, L, C, voltage-source, or current-source dipoles exposed as + current outputs. Positive current follows each netlist ``n1 -> n2`` + declaration. + initial_value: + Initial differential state. Defaults to zero for every state. + + Notes + ----- + Output ports are ordered with all requested voltages first, followed by all + requested currents. Their labels are ``V(node)`` and ``I(dipole)``. + Input ports retain the voltage-source-then-current-source ordering of the + netlist model. The underlying :class:`CircuitModel` is available as + :attr:`model`. + + Examples + -------- + >>> block = NetlistStateSpace( + ... "filter.net", + ... output_voltages=["n1"], + ... output_currents=["Rload"], + ... ) + """ + + def __init__( + self, + netlist: str | Path, + output_voltages: list[str] | None = None, + output_currents: list[str] | None = None, + initial_value: np.ndarray | None = None, + ): + if isinstance(netlist, Path): + elements = parse_netlist_file(netlist) + elif isinstance(netlist, str): + candidate = Path(netlist) + try: + is_file = candidate.is_file() + except OSError: + is_file = False + elements = parse_netlist_file(candidate) if is_file else parse_netlist(netlist) + else: + raise TypeError("netlist must be a string or pathlib.Path") + + self.model = CircuitModel(elements) + for node_name in output_voltages or []: + self.model.add_node_voltage_output(node_name) + for dipole_name in output_currents or []: + self.model.add_dipole_current_output(dipole_name) + + A, B, C, D = self.model.get_system() + super().__init__( + A=A, + B=B, + C=C, + D=D, + initial_value=initial_value, + state_labels=list(self.model.state_labels), + input_labels=list(self.model.input_names), + output_labels=list(self.model.output_labels), + ) diff --git a/tests/filter.net b/tests/filter.net new file mode 100644 index 0000000..a7ce135 --- /dev/null +++ b/tests/filter.net @@ -0,0 +1,16 @@ +*validation netlist for pathsim +.param FREQ = 9k +.param DUTY = 0.5 +.param V_LOW = 0 +.param V_HIGH = 1 +.param T_PER = {1/FREQ} +.param T_ON = {T_PER * DUTY} +.param T_RISE = 10n +.param T_FALL = 10n + +V1 in 0 PULSE({V_LOW} {V_HIGH} 0 {T_RISE} {T_FALL} {T_ON} {T_PER}) +Lfilter in n1 0.5e-3 +Cfilter n1 0 50e-6 +Rload n1 0 50 +.tran 0.00001 0.2 0 0.00001 +.end diff --git a/tests/filter_without_load.net b/tests/filter_without_load.net new file mode 100644 index 0000000..19c47ac --- /dev/null +++ b/tests/filter_without_load.net @@ -0,0 +1,16 @@ +*validation netlist for pathsim +.param FREQ = 9k +.param DUTY = 0.5 +.param V_LOW = 0 +.param V_HIGH = 1 +.param T_PER = {1/FREQ} +.param T_ON = {T_PER * DUTY} +.param T_RISE = 10n +.param T_FALL = 10n + +V1 in 0 PULSE({V_LOW} {V_HIGH} 0 {T_RISE} {T_FALL} {T_ON} {T_PER}) +Bload n1 0 I=V(n1)/50 +Lfilter in n1 0.5e-3 +Cfilter n1 0 50e-6 +.tran 0.00001 0.2 0 0.00001 +.end diff --git a/tests/test_netlist_parser.py b/tests/test_netlist_parser.py new file mode 100644 index 0000000..4faa14c --- /dev/null +++ b/tests/test_netlist_parser.py @@ -0,0 +1,143 @@ +######################################################################################## +## +## TESTS FOR +## 'netlist_to_statespace.py' +## +######################################################################################## + +# IMPORTS ============================================================================== +import unittest +from pathlib import Path +from unittest.mock import patch +import numpy as np + +import pathsim_rf.netlist_to_statespace +from pathsim_rf.netlist_to_statespace import CircuitModel, parse_netlist, parse_value_with_units, parse_netlist_file + +# TESTS ================================================================================ + +class TestNetlistParser(unittest.TestCase): + """Test the netlist parser functions and CircuitModel used in NetlistStateSpace block.""" + + def test_parse_value_with_si_suffixes(self): + """Test conversion between spice si suffixes and power of ten.""" + cases = { + "10f": 10e-15, + "10p": 10e-12, + "10n": 10e-9, + "10u": 10e-6, + "10µ": 10e-6, + "10m": 10e-3, + "10k": 10e3, + "10K": 10e3, + "10meg": 10e6, + "10Meg": 10e6, + "10M": 10e6, + "10g": 10e9, + "10G": 10e9, + "10t": 10e12, + "10T": 10e12, + "100.e-3": 100e-3, + "2.2kOhm": 2200.0, + "3mH": 3e-3, + "4uF": 4e-6, + } + for raw, expected in cases.items(): + self.assertAlmostEqual(parse_value_with_units(raw), expected) + + def test_parse_param_reference(self): + """Test .param in spice list are taken into account""" + elements = parse_netlist( + """ + .param RMAIN = 12k + R1 n1 0 {RMAIN} + L1 n1 0 1mH + """ + ) + self.assertEqual(len(elements), 2) + self.assertAlmostEqual(elements[0].value, 12000.0) + self.assertAlmostEqual(elements[1].value, 1e-3) + + def test_parse_current_source_waveform_keeps_placeholder(self): + """Test if current sources with args are correctly parsed""" + elements = parse_netlist('I1 0 N3 PWL file="SIGNAL.TXT"') + self.assertEqual(len(elements), 1) + self.assertEqual(elements[0].kind, "I") + self.assertIsNone(elements[0].value) + + def test_parse_behavioral_source_placeholder(self): + """Test behavioral sources with are correctly parsed + and segragated between current and voltage sources""" + elements = parse_netlist("B1 0 N3 I=10*(exp(-t)-exp(-2*t)) \n B2 0 N4 V=10*(exp(-t)-exp(-2*t))") + self.assertEqual(len(elements), 2) + self.assertEqual(elements[0].kind, "I") + self.assertIsNone(elements[0].value) + self.assertEqual(elements[1].kind, "V") + self.assertIsNone(elements[0].value) + + def test_unsupported_element_raises(self): + """Test unsupported netlist elements raising an error""" + with self.assertRaises(ValueError): + parse_netlist("X1 n1 0 some_subckt") + + def test_get_system_without_outputs_has_zero_rows(self): + """Test system with zero output correctly formed""" + model = CircuitModel( + parse_netlist( + """ + V1 n1 0 1 + R1 n1 n2 1k + C1 n2 0 1u + """ + ) + ) + A, B, C, D = model.get_system() + self.assertEqual(C.shape, (0, A.shape[0])) + self.assertEqual(D.shape, (0, B.shape[1])) + self.assertEqual(model.output_labels, []) + + def test_outputs_accumulate_for_all_supported_dipoles(self): + """Test output currents can be selected as StateSpace output for all dipole types""" + model = CircuitModel( + parse_netlist( + """ + V1 n1 0 1 + R1 n1 n2 10 + C1 n2 0 1u + L1 n2 0 2m + I1 n2 0 1 + """ + ), + ) + + model.add_node_voltage_output("n2") + for name in ["R1", "L1", "C1", "V1", "I1"]: + model.add_dipole_current_output(name) + A, B, C, D = model.get_system() + + self.assertEqual(C.shape, (6, A.shape[0])) + self.assertEqual(D.shape, (6, B.shape[1])) + self.assertEqual( + model.output_labels, + ["V(n2)", "I(R1)", "I(L1)", "I(C1)", "I(V1)", "I(I1)"], + ) + np.testing.assert_allclose(C[-1], 0.0) + np.testing.assert_allclose(D[-1], [0.0, 1.0]) + + def test_model_without_inputs_raises(self): + """Test error in case of no input present in netlist (in the form of current or voltage source)""" + with self.assertRaisesRegex(ValueError, "at least one independent"): + CircuitModel(parse_netlist("R1 n1 0 1k\nC1 n1 0 1u")) + + def test_unknown_node_and_dipole_raise(self): + """Test error in case of unknown dipole or node selected as output""" + model = CircuitModel(parse_netlist("V1 n1 0 1\nR1 n1 n2 1k\nC1 n2 0 1u")) + with self.assertRaisesRegex(ValueError, "Unknown node"): + model.add_node_voltage_output("missing") + with self.assertRaisesRegex(ValueError, "Unknown dipole"): + model.add_dipole_current_output("R404") + +# RUN TESTS LOCALLY ==================================================================== + +if __name__ == '__main__': + unittest.main(verbosity=2) diff --git a/tests/test_netlist_statespace.py b/tests/test_netlist_statespace.py new file mode 100644 index 0000000..abb6a09 --- /dev/null +++ b/tests/test_netlist_statespace.py @@ -0,0 +1,91 @@ +######################################################################################## +## +## TESTS FOR +## 'netlist_to_statespace.py' +## +######################################################################################## + +# IMPORTS ============================================================================== +import unittest +from pathlib import Path +import numpy as np +import importlib.util +TEST_DIR = Path(__file__).parent + +from pathsim_rf.netlist_to_statespace import CircuitModel, NetlistStateSpace, parse_netlist_file, parse_netlist +from pathsim.blocks.lti import StateSpace + + +# TESTS ================================================================================ + +class TestNetlistStateSpace(unittest.TestCase): + """Test the NetlistStateSpace block""" + + def test_path_string_matches_manual_circuit_model(self): + """Test CircuitModel is correctly used within NetlistStateSpace""" + netlist_path = TEST_DIR / "filter.net" + model = CircuitModel(parse_netlist_file(netlist_path)) + model.add_node_voltage_output("n1") + model.add_dipole_current_output("Rload") + A, B, C, D = model.get_system() + + block = NetlistStateSpace(netlist_path, output_voltages=["n1"], output_currents=["Rload"]) + + self.assertIsInstance(block, StateSpace) + np.testing.assert_allclose(block.A, A) + np.testing.assert_allclose(block.B, B) + np.testing.assert_allclose(block.C, C) + np.testing.assert_allclose(block.D, D) + self.assertEqual(block.state_labels, model.state_labels) + self.assertEqual(block.input_labels, ["V1"]) + self.assertEqual(block.output_labels, ["V(n1)", "I(Rload)"]) + + def test_path_object_preserves_feedback_source_input(self): + """Test current sources feedback used to chain models""" + block = NetlistStateSpace(TEST_DIR / "filter_without_load.net", + output_voltages=["n1"], + output_currents=["Lfilter"]) + + self.assertEqual(block.input_labels, ["V1", "Bload"]) + self.assertEqual(block.output_labels, ["V(n1)", "I(Lfilter)"]) + self.assertEqual(len(block.inputs), 2) + self.assertEqual(len(block.outputs), 2) + + def test_inline_netlist_string(self): + """Test direct netlist use to class for instance generation""" + block = NetlistStateSpace( + """ + V1 in 0 1 + R1 in out 10 + C1 out 0 1u + """, + output_voltages=["out"], + output_currents=["R1", "V1"]) + + self.assertEqual(block.input_labels, ["V1"]) + self.assertEqual(block.output_labels, ["V(out)", "I(R1)", "I(V1)"]) + self.assertEqual(block.C.shape[0], 3) + self.assertEqual(block.D.shape[0], 3) + + def test_explicit_missing_path_raises(self): + """Test missing .net file provided""" + with self.assertRaises(FileNotFoundError): + NetlistStateSpace(TEST_DIR / "missing.net") + + def test_path_str_import(self): + """Test str netlist_path arg allows file to be imported""" + netlist_path = TEST_DIR / "filter.net" + NetlistStateSpace(str(netlist_path)) + + def test_invalid_output_names_raise(self): + """Test invalid output names provided for state space generated""" + netlist_path = TEST_DIR / "filter.net" + with self.assertRaisesRegex(ValueError, "Unknown node"): + NetlistStateSpace(netlist_path, output_voltages=["missing"]) + with self.assertRaisesRegex(ValueError, "Unknown dipole"): + NetlistStateSpace(netlist_path, output_currents=["R404"]) + + +# RUN TESTS LOCALLY ==================================================================== +if __name__ == '__main__': + unittest.main(verbosity=2)