Control Systems
A transfer function in Leonardo is an ordinary expression — a Ratio of polynomials in a variable you name. There is no TransferFunction type, no model object, and nothing to construct: 1/(s+1) is one, and so is anything the rest of the library produces that happens to be rational.
That is a deliberate decision, and it has a visible consequence in every signature on this page: the frequency variable is always an explicit argument, because a plain Ratio carries no type to dispatch on. You write poles(g, s), never g.poles.
import it.grypho.scala.leonardo.core.*
import it.grypho.scala.leonardo.control.*
import it.grypho.scala.leonardo.parser.Parser
val env = new Environment()
val s = _Variable("s")
val t = _Variable("t")
val z = _Variable("z")
def tf(src: String): _Expression = Parser.parse(src).get
Why there is no carrier
The argument for a dedicated type is real: it would let the frequency variable be implicit and would give G * H a meaning. The argument against it won, on three counts.
A carrier would be a wall. simplify, derive, integrate, substitute, the exact tier and the parser all operate on expressions; a transfer function that is not one would need each of them taught about it, or would be cut off from all of them. As an expression, a transfer function is differentiable, simplifiable, exactly representable and printable for free.
The one thing a carrier genuinely buys — a result that is a single tidy fraction rather than a nested one — is supplied here instead: every interconnection returns N(s)/D(s) in lowest terms.
And the cost of the decision is paid in documentation, which is what this page is.
Interconnection
Three connections, each returning one normalised rational.
series(tf("1/(s+1)"), tf("2/(s+3)"), s).toString
// res0: String = "(2.0 / ((3.0 + (4.0 * s)) + (s ^ 2.0)))"
parallel(tf("1/(s+1)"), tf("1/(s+2)"), s).toString
// res1: String = "((3.0 + (2.0 * s)) / ((2.0 + (3.0 * s)) + (s ^ 2.0)))"
feedback is negative feedback, G/(1 + G·H):
feedback(tf("1/(s*(s+2))"), _Number(1), s).toString
// res2: String = "(1.0 / ((1.0 + (2.0 * s)) + (s ^ 2.0)))"
The sign convention is stated because it cannot be inferred. Negative feedback is the control convention and what a unity-feedback loop means by default, but a substantial part of the literature writes the positive form G/(1 − G·H). For positive feedback, negate h. This is the same care the spherical polar angle convention needed: choosing silently would make every closed-loop result wrong for half its readers.
Notice what the closed loop above returned: 1/(s² + 2s + 1), not (s² + 2s)/((s² + 2s)(s² + 2s + 1)). The common factor is divided out, because otherwise the poles would come back as {0, −2, −1, −1} — a confident answer to a different question.
Poles, zeros and stability
poles(tf("(s+2)/((s+1)*(s+3))"), s).map(_.mkString(", "))
// res3: Option[String] = Some(value = "-1.0, -3.0")
zeros(tf("(s+2)/((s+1)*(s+3))"), s).map(_.mkString(", "))
// res4: Option[String] = Some(value = "-2.0")
dcgain(tf("5/(s+2)"), s)
// res5: Option[Double] = Some(value = 2.5)
Poles are _Values, not Doubles, because a complex pole is the interesting case — an oscillatory mode is a conjugate pair, and reducing the result to reals would silently drop exactly the systems control theory is about.
poles(tf("1/(s^2 + 1)"), s).map(_.mkString(", "))
// res6: Option[String] = Some(value = "i, -i")
Stability is decided from the poles, strictly:
isStable(tf("1/(s^2 + 3*s + 2)"), s)
// res7: Option[Boolean] = Some(value = true)
isStable(tf("1/(s - 1)"), s)
// res8: Option[Boolean] = Some(value = false)
isStable(tf("1/(s^2 + 1)"), s)
// res9: Option[Boolean] = Some(value = false)
The last one is the marginal case, and it is false. An oscillator that never decays is not a stable system. A Routh array is available as routhTable(g, s) for inspecting a margin by hand — from Scala only, since it has no grammar production and routh is consequently an ordinary variable name, not a reserved word — but isStable does not consult it: Routh’s degenerate cases (a zero in the first column, an identically zero row) each need their own repair, and a mishandled one yields a wrong verdict rather than a refusal. Root location has neither failure mode.
Time response
stepResponse(tf("1/(s+1)"), s, t).toString
// res10: String = "(1.0 + (-1.0 * exp((-1.0 * t))))"
impulseResponse(tf("2/(s+3)"), s, t).toString
// res11: String = "(2.0 * exp((-3.0 * t)))"
These are the inverse Laplace transform of G(s)/s and of G(s) — almost entirely a re-spelling of machinery that already existed. The response variable t is free in the result, so it can be evaluated, plotted with sample, or differentiated like any other expression.
Frequency response
bode(tf("1/(s+1)"), s, 1.0)
// res12: Option[Tuple2[Double, Double]] = Some(
// value = (0.7071067811865476, -0.7853981633974483)
// )
nyquist(tf("1/(s+1)"), s, 1.0)
// res13: Option[Tuple2[Double, Double]] = Some(
// value = (0.5000000000000001, -0.5)
// )
(magnitude, phase-in-radians) — the data of a Bode plot — and (real, imaginary) — the data of a Nyquist plot — one computation presented two ways. At the corner frequency of a first-order lag the magnitude is 1/√2 (−3 dB) and the phase −45°, as it should be. The whole of the frequency response is a single substitution s → iω riding the existing complex closure; no new arithmetic was written for it.
An integrator has no finite response at ω = 0, and that is reported as absence rather than a fabricated infinity:
bode(tf("1/s"), s, 0.0)
// res14: Option[Tuple2[Double, Double]] = None
Sweeping a band
bode and nyquist answer for one frequency. A diagram needs a vector, and it needs two things that a plain sample cannot give it:
frequencyResponse(tf("1/(s+1)"), s, 0.1, 10.0, 5, new Environment())
// res15: Vector[Tuple3[Double, Double, Double]] = Vector(
// (0.1, -0.04321373782642559, -5.710593137499643),
// (0.31622776601683805, -0.41392685158225007, -17.548400613792307),
// (1.0000000000000004, -3.0102999566398143, -45.00000000000001),
// (3.162277660168381, -10.413926851582254, -72.45159938620772),
// (10.0, -20.043213737826427, -84.28940686250037)
// )
(ω, magnitude in dB, phase in degrees) — dB and degrees because that is what a Bode plot is; the raw pair stays available from bode itself.
The grid is geometric, not linear. A frequency response is read across decades, so a linear grid of 200 points over 0.01 .. 100 would put 199 of them in the last decade and none near a corner at 0.1 — the interesting part of the curve is exactly the part it fails to resolve.
The phase is unwrapped, and this is the substance. atan2 has principal value (−π, π], so a swept curve jumps a full turn wherever it crosses the branch cut. That jump is an artefact of the arctangent, not of the plant. A third-order lag tends to −270°:
frequencyResponse(tf("1/(s+1)^3"), s, 1.0, 100.0, 3, new Environment()).map(_._3)
// res16: Vector[Double] = Vector(
// -135.0,
// -252.86822058750104,
// -268.2811839069496
// )
A raw sweep would report the last of those as +90. One caveat worth stating: unwrapping cannot tell a genuine half-turn step from a grid too coarse to resolve a fast one, so a sparse sweep across a lightly damped resonance can unwrap the wrong way — the remedy is points.
nyquistSweep reads the same grid and returns (real, imaginary), the one sweep presented two ways that bode and nyquist already are.
Both take an Environment, so a plant carrying a bound parameter folds — which is what lets the browser REPL draw a diagram from a session where you wrote K := 10 first. There, bode G s 0.01 100 and nyquist G s 0.01 100 draw them.
State space
A state-space model is a 1×4 row of matrices — the same shape lu, qr and eig return — so at(m, 1, k) indexes it and a session :save round-trips it with no new machinery. Controllability asks whether the input can steer every state; observability whether the output can distinguish them:
controllable(tf("[[0, 1], [-2, -3]]"), tf("[[0], [1]]"))
// res17: Option[Boolean] = Some(value = true)
controllable(tf("[[1, 0], [0, 2]]"), tf("[[1], [0]]"))
// res18: Option[Boolean] = Some(value = false)
Observability is defined as controllability of the dual pair (Aᵀ, Cᵀ) rather than restated, so the two can never disagree.
Both answers are independent of the units the model is written in, which is less obvious than it sounds. Rank is decided by QR of the transposed controllability matrix, not by det(M·Mᵀ) against a threshold: a determinant scales like ‖M‖^(2n), so with the older test the same plant with B in millivolts rather than volts came back uncontrollable. Scaling B cannot change which states the input can reach, and the tests pin that it does not.
Discrete time
Discretisation always names its method — the zero-order hold Zoh or the bilinear (Tustin) transform Tustin. The same plant discretised the two ways has different discrete poles, so a silent default would make two correct-looking answers disagree with no way to see why.
c2d(tf("1/(s+1)"), s, z, 0.1, Tustin).map(_.toString)
// res19: Option[String] = Some(value = "((1.0 + z) / (-19.0 + (21.0 * z)))")
Tustin is exactly invertible, which is what makes the round trip return to where it started. The result below is 2/(s+3) with a factor of -0.2 left in both numerator and denominator — Leonardo does not cancel a common numeric factor out of a fraction, so read the value, not the spelling:
val back = c2d(tf("2/(s+3)"), s, z, 0.05, Tustin).flatMap(d2c(_, z, s, 0.05, Tustin))
back.map(_.toString)
// res20: Option[String] = Some(value = "(-0.4 / (-0.6 + (-0.2 * s)))")
back.map(_.eval(env.withBinding("s", _Number(1.0)))) // 2/(1+3) = 0.5
// res21: Option[Either[_Expression, _Value]] = Some(
// value = Right(value = _Number(d = 0.5000000000000003))
// )
Discrete stability is the unit circle, not the left half-plane — applying the continuous rule to a z-domain function inverts the answer rather than degrading it:
isStableDiscrete(tf("1/(z - 0.5)"), z)
// res22: Option[Boolean] = Some(value = true)
isStableDiscrete(tf("1/(z - 1.5)"), z)
// res23: Option[Boolean] = Some(value = false)
Exact state-space discretisation uses the matrix exponential of the block matrix [[A, B], [0, 0]]·Ts, not B_d = A⁻¹(A_d − I)B. That is a correctness matter: the block form is defined for a singular A, and any system with an integrator has one.
c2dExact(tf("[[0, 1], [0, 0]]"), tf("[[0], [1]]"), 0.5)
// res24: Option[Tuple2[_MatrixValue, _MatrixValue]] = Some(
// value = ([[1.0, 0.5], [0.0, 1.0]], [[0.125], [0.5]])
// )
What it declines
The refusals are as much of the design as the features.
A dead-time term is not rational, so it has no poles — and it is not silently Padé-approximated. pade exists, and reaching for it is the user’s decision:
poles(tf("exp(-2*s)/(s+1)"), s).map(_.mkString(", "))
// res25: Option[String] = None
An improper G has a step response containing an impulse at t = 0 that the inverse-transform tier does not represent, so it returns unevaluated rather than handing back only the smooth part:
stepResponse(tf("s^2/(s+1)"), s, t).toString
// res26: String = "step(((s ^ 2.0) / (s + 1.0)), s, t)"
A coefficient whose sign cannot be determined leaves stability undecided, rather than guessed:
isStable(tf("1/(s^2 + a*s + 1)"), s)
// res27: Option[Boolean] = None
Also absent, deliberately: d2c by zero-order hold (it needs a matrix logarithm), root-locus and margin plots (the numbers are here; plotting is not this library’s job), and a discrete Kronecker delta, which core has no node for.
From the REPL
> feedback(1/(s*(s+2)), 1, s)
> step(1/(s+1), s, t)
> impulse(1/(s+1), s, t)
Note that step is arity-overloaded: step(x) is the Heaviside unit step, while step(G, s, t) is the step response of a system.