Integration with python-control

[milanofthe]

@murrayrm, @RemDelaporteMathurin lets move the discussion about collaboration/integration here.

Richard and I had a bit of an email exchange (I was reaching out to the python-control team), I will paste the core here so we can continue the discussion in the open.

His ideas were:

• "Exporting" python-control models to pathsim: given an InterconnectedSystem (our version of a block diagram), we could chase through the structure and call pathsim to create the equivalent pathsim blocks (lti, ode, etc) and connections.
• Linearizing pathsim models to be able to carry out linear controller designs (LQR, hinf, etc). We would then want to export the controller (either static gains or a dynamic compensator) as a pathsim block.
• Wrapping pathsim models with a NonlinearIOSystem to allow them to be used in python-control. This would probably be hard, but would allow the use of the optimal control module and other functions that do more than just simulate the system.

My response to them was:

• The export thing makes a lot fof sense for cross compatibility. I do have some dedicated control blocks in the ctrl module but having this done through python-control or as an addition probably has some value.
• I was talking to some control folks about this some months ago. PathSim does support block level linearization, but this would require end-to-end linearization of whole systems ore subsystems to get linear state space matrices. I would need to think about a good way of doing that.
• I agree, doable in principle. The difficulty primarily depends on what interface to the system python control needs and in what form (time/s domain). I am not super familiar there. Currently I have a contributor working on a parameter fitting and optimization module for PathSim. I will talk to him, maybe there is some overlap.

And his gain:

Do you happen to know anyone who is contributing to pathsim that is familiar with python-control? If so, I think they would be in a good position to add a to_pathsim method to the various InputOutputSystem subclasses that would do the trick. We already have something like for generating Pandas dataframes when Pandas is present.
In terms of the linearization, you could possibly get away with a numerical solution, similar to what we do in the NonlinearIOSystem.linearize method. This would only work around an operating point and would have to have constant values for any discrete/switching variables, but could be useful.

[miguelpi314]

Hi @milanofthe (and @murrayrm if still following this thread),

I'd like to pick up the "end-to-end linearization of whole systems" gap mentioned above. Looking at the current code, the building blocks are already there but not connected:

• Simulation.steadystate() finds a pseudo-equilibrium (forces ẋ=0 at fixed inputs), but isn't a general trim — no way to fix target outputs and solve for the inputs/states that produce them.
• Simulation.linearize() / Subsystem.linearize() already linearize every block's internal Operator/DynamicOperator individually, but nothing walks the connection graph afterward to assemble a single global (A, B, C, D) for the interconnected system.
• StateSpace in blocks/lti.py is already the natural container for that result.

I'd like to propose closing that gap with:

1. Simulation.trim(spec) — generalize steadystate() to accept a mix of fixed states/inputs/outputs and solve for the rest, reusing the existing Newton/Anderson optimizer in optim/.
2. Simulation.linearize_system() — walk blocks + connections, collect the per-block Jacobians already computed by block.linearize(), and assemble them into a single (A, B, C, D) respecting external ports (handling algebraic loops the same way the solver already does internally), returned as a StateSpace block or raw matrices.
3. A worked example (PID + nonlinear plant) showing trim → linearize → classical analysis (stability margins, root locus) using plain matplotlib, no external dependency required for the MVP.

I see this as the foundation the other two ideas in this thread depend on: once there's a reliable global (A,B,C,D), exporting it into python-control (or the reverse, a to_pathsim() importer for InterconnectedSystem) becomes much more tractable, and so does eventually wrapping a pathsim model as a NonlinearIOSystem. I'd rather get the linearization/trim piece solid first than try to tackle all three at once.

Happy to scope this as a PR against pathsim. Let me know what do you think before I start.

[milanofthe]

Hi @miguelpi314 awesome that you are willing to take this on!

[milanofthe]

Picking this thread back up, because @miguelpi314's PR #237 made me work out what
the right shape actually is. Writing the design down here so it's in the open
before anyone builds against it.

Short version: linearization in PathSim is becoming a query that returns
matrices, not a mode the system is switched into. That's a change to something
that has been in the core for a while, so here's the reasoning.

Why the current design doesn't extend

Today Operator.linearize(x) / DynamicOperator.linearize(x, u, t) compute the
Jacobian and store it, after which __call__ returns the first-order Taylor
surrogate until reset(). Block.linearize() / Simulation.linearize() are the
switches on top. The Jacobian is a by-product living inside the operator, not
something you can ask for.

That works fine for what it was built for — flip the system to linear, simulate,
flip back. It does not extend to assembling a global (A, B, C, D), because
assembly needs to read every block's Jacobian without changing how the system
behaves afterwards. Building the read on top of the switch means the simulation
silently keeps integrating the surrogate once you've asked for a model. That is
what happens in #237, and it's not really a bug in the PR so much as the
consequence of the layer it had to build on.

Where fastsim already landed

fastsim hit the same question earlier and answered it the
other way. From src/blocks/operator.rs:

the Jacobian, derived by forward/reverse AD from the graph itself
(sparse-aware, gated, width-cached) — NOT a hand-written callback, so it can
never desync from the function … mirroring pathsim's Operator(func, jac)
abstraction (but the AD is over our SSA graph, and the Jacobian is recomputed
per Newton step rather than linearized-and-reused).

The Jacobian there is Operator::jacobian_wrt_state(x, u, t, mem) — a pure
query, derived from the SSA graph, never a stored mode. Simulation.linearize()
on that side is deliberately NotImplementedError.

I want the two projects to agree here, so PathSim moves to the same shape.

The design

Three levels, to_statespace() at each:

1. Operator — already done. Operator.jac(x) and
DynamicOperator.jac_x(x, u, t) / .jac_u(x, u, t) are already pure,
side-effect-free queries. blocks/dae.py uses op_dyn.jac_x() in its solve path
today. Nothing to build; the layer just needs to be the one that gets used.

2. Block — the actual gap. Block.to_statespace(t) returns the block's local
linear model as a StateSpace. The base implementation reads the operator
queries. Blocks without an operator override it with exact values rather than
finite differences — an Integrator is A=0, B=I, C=I, D=0 and there is no
reason to difference that numerically.

Blocks that carry signal but have no valid linear model — Comparator, Relay,
Switch, RNG, Counter, the discrete and logic blocks — raise from their
linearize(). A blanket numerical fallback produces a "Jacobian" for a
comparator at its switching point, which is worse than an error. This mirrors how
fastsim separates traceable operators from genuinely opaque ones.

Note what that implies: Simulation.linearize() on a system containing such a
block starts raising where it previously succeeded silently. That is intended. A
"linearized" system with a relay in it was never a linear system, and PathSim is
pre-1.0 — better to fix the semantics now than to carry them past a stable
release. Pure sinks like Scope are unaffected: they have no outputs, so their
linear model is legitimately empty rather than undefined.

Block.linearize(t) otherwise keeps doing what it does today and additionally
returns the matrices it computes anyway, so the switch becomes a thin wrapper
over the query and existing code keeps working.

3. Simulation / Subsystem. Simulation.to_statespace(inputs, outputs, t)
assembles the global model. I want this done in the interconnection form —
block-diagonal (A, B, C, D) plus an interconnection matrix L, solved as

v = (I - L·D)⁻¹ (L·C·x + M·u)

rather than by substituting along a topological order. Two reasons: it handles
systems where an algebraic loop survives the input break, which PathSim's own
_update() already solves and a DAG walk has to reject; and well-posedness
becomes the invertibility of (I - L·D) instead of a special case. It is also
considerably less code.

Since every level returns a StateSpace block, Subsystem.to_statespace() falls
out for free, and hierarchical models linearize hierarchically.

What this gives the python-control side

@murrayrm — mapping this back to your three ideas:

• Linearizing PathSim models (Packaging improvements #2 on your list): this is it. The output is a
  StateSpace with A/B/C/D plus state/input/output labels, which construct a
  control.StateSpace directly. Add trim() and linearize_system() for system-level operating points #237 verified that part against real
  python-control — poles, DC gain and step response all agree — and that
  verification stands regardless of how the matrices get assembled. Exporting a
  designed controller back is then just constructing a StateSpace block, which
  already works.
• Export from python-control (PID with anti-windup implementation #1): independent of this, and probably the
  easier one. An InterconnectedSystem walker emitting PathSim blocks and
  connections doesn't need any of the above.
• Wrapping as NonlinearIOSystem (Warn if block is in .connections but not in simulation.blocks #3): still the hard one, and still open. I
  don't want to design it as part of this.

Trim

The other half of #237 is a trim() that solves backwards for the input
producing a desired output. I want the capability — steadystate() genuinely
can't do it — but as an extension of steadystate() rather than a second solver
stack next to it. steadystate() works by switching every block to the
SteadyState solver and reusing the existing _solve(); trim is that same solve
with the unknown vector extended by the free parameters and the residual extended
by the output targets.

That also answers the question of whether trim is provisional: it stays PathSim's
own, because it's an extension of machinery PathSim already has rather than a
thin thing that a find_operating_point() bridge would supersede.

Status

Level 1 exists. Levels 2 and 3 are open — level 2 is self-contained and a good
first piece if anyone wants it. @miguelpi314, the offer from the PR thread stands.