Python API

ODE and DAE API

Symbolic ordinary and differential-algebraic equation systems.

ODE

class

class ODE(state_vars, rhs, time_var)

#source

Symbolic system of first-order ordinary differential equations.

Represents the system d(state_vars)/dt = rhs driven by time_var.

Parameters

state_varslist[Expr]: Symbolic state variables (e.g. [x, y]).
rhslist[Expr]: Right-hand side expressions, one per state variable.
time_varExpr: The independent time variable.
Example:: p = alkahest.ExprPool() t, x = p.symbol("t"), p.symbol("x") ode = alkahest.ODE([x], [-x], t) # dx/dt = -x

Properties

order

#source

Number of state variables (the order of the system).

Methods

new(state_vars, rhs, time_var)

staticmethod

#source

Build an explicit first-order system d(state_vars)/dt = rhs.

Parameters

state_varslist[Expr]: Symbolic state variables.
rhslist[Expr]: Right-hand sides, one per state variable and in the same order.
time_varExpr: The independent variable.

All three arguments are positional; there are no keyword forms. ODE(state_vars, rhs, time_var) is accepted as a synonym. For an implicit or constrained system use DAE instead, and for a scalar higher-order equation use lower_to_first_order().

Example:

Python

# Harmonic oscillator as a first-order system: x' = v, v' = -x
ode = alkahest.ODE.new([x, v], [v, -x], t)

from_equations(equations, functions, time_var)

staticmethod

#source

Build the first-order system equivalent to equations written with applied functions.

Parameters

equationslist[Expr]: One expression per unknown, each read as = 0, written with y(t) and diff(y(t), t) (any order). They must be affine in the highest derivative of each unknown, with an invertible coefficient matrix.
functionslist[Expr]: The unknowns, [y(t), z(t), ...] (undefined functions of t).
time_varExpr: The independent variable t.

The state of a function with highest derivative of order m is [y, y_dot, ...] (m symbols named like functions_to_symbols() names them). Refusals: E-DIFF-005/007 for an unknown that is not an undefined function of t; E-ODE-030 when the system does not reduce (wrong equation count, an unknown that is never differentiated, a nonlinear or singular leading part, i.e. a DAE). A leading coefficient that needs != 0 is refused here; pass those equations to dsolve_system(), which reports the condition.

with_ic(var, value)

#source

is_autonomous()

#source

state_vars()

#source

rhs()

#source

derivatives()

#source

simplify_rhs()

#source

experimental.dsolve_system(ode_or_equations, assumptions=None, *, functions=None, t=None, real_form=True) Closed-form solution of y' = A y + f(t). Accepts an ODE or a list of equations in y(t) form (then functions and t are required, and the result carries "functions": {y(t): expr}). With real_form the closed form of a real system is written with cos/sin/cosh/sinh wherever the imaginary part is proved zero (see to_real_form()).

function

lower_to_first_order(var, rhs, order, time_var)

#source

alkahest.lower_to_first_order(var, rhs, order, time_var) — lower a scalar higher-order ODE to an equivalent first-order system.

Takes the pieces of the equation, not an ODE object: the scalar unknown, the right-hand side, the order of the derivative on the left, and the independent variable. All four are required and positional.

Parameters

varExpr: The scalar unknown, e.g. x.
rhsExpr: Right-hand side of d^order(var)/dt^order = rhs.
orderint: Order of the equation (2 for x'').
time_varExpr: The independent variable.

Returns

ODE: A first-order system whose state is [var, var', …, var^(order-1)].
Example:: ode = alkahest.lower_to_first_order(x, -4 * x, 2, t) ode.order # 2

DAE

class

class DAE

#source

Symbolic differential-algebraic equation system.

A DAE is a system of implicit equations mixing differential and algebraic constraints: F(t, variables, derivatives) = 0.

Parameters (via :meth:`DAE.new`)

equations : list[Expr] Implicit equations (each equals zero). variables : list[Expr] Algebraic and state variables. derivatives : list[Expr] Derivative expressions corresponding to variables. time_var : Expr The independent time variable.

Use pantelides() to reduce the differential index before simulation.

Attributes

n_equationsint
n_variablesint
time_varExpr
indexint | None: Differentiation index, set by pantelides() on the DAE it returns; None on a DAE that has not been index-reduced.
Example:: p = alkahest.ExprPool() t, x, dx = p.symbol("t"), p.symbol("x"), p.symbol("dx/dt") dae = alkahest.DAE.new([dx - x], [x], [dx], t) # x' = x dae.equations() # [-x + dx/dt]

Properties

n_equations

#source

Number of equations in the system.

n_variables

#source

Number of dependent variables in the system.

time_var

#source

The independent (time) variable.

index

#source

Differentiation index, or None if it has not been computed.

pantelides() sets it on the DAE it returns: 0 for an explicit ODE, 1 for a semi-explicit index-1 system, 3 for the Cartesian pendulum — the number of differentiations needed to determine every derivative, including one for the algebraic variables.

Methods

new(equations, variables, derivatives, time_var)

staticmethod

#source

Build a DAE from implicit equations.

Parameters

equationslist[Expr]: Implicit equations, each meaning g = 0. Write x' = f as dx - f, using a separate symbol for the derivative.
variableslist[Expr]: Dependent variables (state and algebraic).
derivativeslist[Expr]: Symbol standing for the time derivative of variables[i], e.g. pool.symbol("dx/dt"). Alkahest does not parse the name; the positional pairing is what makes it a derivative. It may be shorter than variables: the trailing variables are then algebraic (their derivative, if ever needed, is named d{name}/dt). A variable listed as another's derivative is a jet of that unknown, not a new one.
time_varExpr: The independent variable.

Raises

DaeError: E-DAE-004 when derivatives is longer than variables, an entry is not a symbol, a symbol is repeated, or time_var is listed among the unknowns.
Example:: x' = u, y' = w, u' = -lamx, w' = -lamy - 1, x2 + y2 = 1

dae = alkahest.DAE.new( [dx - u, dy - w, du + lam * x, dw + lam * y + one, x2 + y2 - one], [x, y, u, w, lam], [dx, dy, du, dw], t, )

equations()

#source

The implicit equations, each meaning g = 0.

After pantelides() or rosenfeld_groebner() this is how you see which equations were added by differentiation.

variables()

#source

The dependent variables, parallel to derivatives().

derivatives()

#source

The derivative symbols: derivatives()[i] is d(variables()[i])/dt.

Index reduction appends higher jets here (d2x/dt2, …), so this is how the extra variables of a prolonged system get names.

function

pantelides(dae)

#source

alkahest.pantelides(dae) — structural index reduction of a DAE.

Matches equations to unknowns on their highest derivatives (Pryce's Σ-method, equivalent to Pantelides' algorithm), differentiates each equation as often as that matching requires, and returns the reduced DAE — the original equations followed by their derivatives — not a report object. Read them with DAE.equations() and the new jets with DAE.derivatives().

DAE.index on the result is the differentiation index: 0 for an explicit ODE, 1 for a semi-explicit index-1 system, 2 for x' = y, x = sin t, 3 for the Cartesian pendulum.

Before returning, the Jacobian of the reduced system with respect to the highest derivatives is checked to be non-singular at generic sample points (exact rational arithmetic); a structural answer the coefficients contradict is refused, never returned.

Raises

DaeError: (a ValueError subclass) E-DAE-002 when an equation would need more than 10 differentiations — try dae_index_reduce(), which falls back to rosenfeld_groebner(); E-DAE-003 when there is not one equation per unknown function or no matching exists; E-DAE-004 for a malformed DAE; E-DAE-005 when the reduced Jacobian is singular.
Example:: reduced = alkahest.pantelides(dae) reduced.index # 3 for the Cartesian pendulum reduced.equations() # original equations plus the differentiated ones

Sensitivity analysis

function

sensitivity_system(ode, params)

#source

alkahest.sensitivity_system(ode, params) — build the forward sensitivity ODE.

function

adjoint_system(ode, objective_grad)

#source

alkahest.adjoint_system(ode, objective_grad) — build the adjoint ODE.

class

class SensitivitySystem

#source

The extended ODE with sensitivity equations attached.

Properties

extended_ode

#source

original_dim

#source

n_params

#source

Acausal modeling

class

class AcausalSystem(pool)

#source

Acausal component-based modelling system.

An AcausalSystem aggregates components connected through Port objects (potential/flow pairs). Add components with add_component(), wire ports together with connect(), and call flatten() to convert the component network into an equivalent DAE suitable for simulation or index reduction.

Example (RC circuit):

Python

p = alkahest.ExprPool()
t = p.symbol("t")

src = alkahest.voltage_source("V1", p.symbol("Vs"))
res = alkahest.resistor("R1", p.symbol("R"))
cap = alkahest.capacitor("C1", p.symbol("C"))

sys = alkahest.AcausalSystem(p)
sys.add_component(src["component"])
sys.add_component(res["component"])
sys.add_component(cap["component"])

sys.connect(src["component"].port("V1.p"), res["component"].port("R1.p"))
sys.connect(res["component"].port("R1.n"), cap["component"].port("C1.p"))
sys.connect(cap["component"].port("C1.n"), src["component"].port("V1.n"))

dae = sys.flatten(t)

Methods

add_component(component)

#source

Add a component (e.g. from resistor(), capacitor(), voltage_source()) to the system.

connect(port_a, port_b)

#source

Connect two ports: equates their potentials and balances their flows (a.potential == b.potential, a.flow + b.flow == 0).

flatten(time_var)

#source

Flatten all component and connection equations into a DAE.

class

class Component

#source

A physical component (resistor, capacitor, voltage source, …) with named Port connectors and internal constitutive equations.

Construct components via resistor(), capacitor(), or voltage_source(), then register them on an AcausalSystem with AcausalSystem.add_component() and wire them together with AcausalSystem.connect().

Properties

name

#source

n_equations

#source

Number of constitutive equations contributed by this component.

n_ports

#source

Number of external connection ports.

Methods

ports() → list[Port]

#source

All ports, in declaration order.

port(name)

#source

Look up a port by its full name (e.g. "R1.p"), or None if absent.

class

class Port

#source

A connection port on a component (e.g. positive/negative terminals of a resistor), exposing potential and flow Expr attributes.

Properties

potential

#source

flow

#source

name

#source

function

resistor(name, resistance)

#source

alkahest.resistor(name, resistance) — create a resistor component.

Returns a dict {"name", "n_equations", "n_ports", "component"}; the "component" entry is a Component usable with AcausalSystem.add_component and .port(name).

function

capacitor(name, capacitance)

#source

alkahest.capacitor(name, capacitance) — create an ideal capacitor component (C * dv/dt = i).

Returns a dict {"name", "n_equations", "n_ports", "component"}; the "component" entry is a Component usable with AcausalSystem.add_component and .port(name).

function

voltage_source(name, voltage)

#source

alkahest.voltage_source(name, voltage) — create an ideal voltage source component (v_p - v_n = V).

Returns a dict {"name", "n_equations", "n_ports", "component"}; the "component" entry is a Component usable with AcausalSystem.add_component and .port(name).

Hybrid systems

class

class HybridODE

#source

Ordinary differential equation system with discrete events.

A HybridODE wraps a continuous ODE and a list of Event objects. Each event specifies a guard condition and a reset map that is applied when the guard crosses zero.

Construction:

Python

hybrid = alkahest.HybridODE.new(ode)
hybrid = hybrid.add_event(event)

Attributes

n_eventsint: Number of discrete events attached to this system.

Properties

n_events

#source

Number of registered discrete events.

Methods

new(ode)

staticmethod

#source

add_event(event)

#source

guards()

#source

class

class Event

#source

A hybrid system event.

Parameters

condition: Expression; event fires when this crosses zero.
reset: Dict mapping state variables to their post-event values.

Methods

new(name, condition, reset_map)

staticmethod

#source

rising()

#source

falling()

#source

Differential algebra

function

dae_index_reduce(dae, order=None)

#source

No docstring.

class

class DaeIndexReduction

#source

No docstring.

Properties

used_pantelides

#source

used_rosenfeld_groebner

#source

Methods

dae()

#source

Pantelides-reduced DAE if Pantelides succeeded; else Rosenfeld working DAE.

rosenfeld_groebner_result()

#source

function

rosenfeld_groebner(
    dae,
    order=None,
    max_prolong_rounds=None,
    params=None,
    eliminate=None,
    minimal=False
)

#source

alkahest.rosenfeld_groebner(dae, order=None, max_prolong_rounds=None) — Rosenfeld–Gröbner-style differential elimination.

Prolongs the system (differentiates each equation, introducing higher jets as new indeterminates) and computes a Gröbner basis after each round, until differentiating adds nothing new to the ideal or the round budget runs out. The basis is the set of algebraic consequences of the differential system — the input–output relations elimination is after.

Parameters

daeDAE

The system, polynomial in its variables and derivative symbols.

orderstr, optional

Monomial order — "grevlex" (default), "grlex" or "lex". Use "lex" when you want elimination-ordered generators.

max_prolong_roundsint, optional

Prolongation budget (default 8). Nonlinear jets often do not saturate in finitely many algebraic steps, so hitting the budget is normal and sets RosenfeldGroebnerResult.truncated.

paramslist[Expr], optional

Symbols to put in the coefficient field Q(params) rather than in the ring (M9). Without this every free symbol is a ring variable, so a model parameter enlarges the monomial order, the pair schedule and the staircase — which is the difference between eliminating states from Q[states, jets, params] and from Q(params)[states, jets]. With params the return type is ParametricRosenfeldGroebnerResult and the basis is ParametricGroebnerBasis; without it nothing changes.

The parameters must not be DAE variables, derivatives or the time variable.

eliminatelist[Expr], optional

Variables the caller intends to eliminate — typically the unobserved states. Their whole jet chain goes with them: naming x also names dx/dt, d2x/dt2, … as prolongation introduces them. Supplying this makes ParametricRosenfeldGroebnerResult.minimal_prolongation_rounds available and enables the over-prolongation warning. Requires params.

minimalbool, optional

Stop at the first prolongation round whose elimination ideal with respect to eliminate is non-empty, instead of prolonging to the budget. One prolongation too many is expensive out of all proportion — on SIR it is 0.002 s and one 4-term relation against 20.9 s and thirteen generators of up to 233 terms — and the over-supplied result is correct, so nothing else signals it. Requires eliminate.

Not a minimality guarantee: see ParametricRosenfeldGroebnerResult.minimal_prolongation_rounds for the scope, in particular for multi-output models.

Returns

RosenfeldGroebnerResult or ParametricRosenfeldGroebnerResult: Read the relations with result.final_basis().to_exprs().

Warns

UserWarning

When eliminate is given, the elimination ideal was already non-empty at an earlier round, and minimal was not set.

Example:

t, x, dx = p.symbol("t"), p.symbol("x"), p.symbol("dx/dt") dae = alkahest.DAE.new([dx - x], [x], [dx], t) r = alkahest.rosenfeld_groebner(dae, max_prolong_rounds=1) r.consistent # True r.final_basis().to_exprs() # the eliminated relations, as Expr

dae = alkahest.DAE.new([dx - a*x], [x], [dx], t) r = alkahest.rosenfeld_groebner(dae, params=[a], max_prolong_rounds=1) r.final_basis().conditions()

class

class RosenfeldGroebnerResult

#source

V2-13 — result of Rosenfeld–Gröbner-style differential elimination.

Returned by rosenfeld_groebner(). The eliminated relations are in final_basis(); read them with result.final_basis().to_exprs().

Attributes

consistentbool: False iff the unit ideal was reached — the system has no common jet solution over ℚ, i.e. the equations are contradictory.
truncatedbool: True if prolongation stopped at max_prolong_rounds rather than because the differential chain saturated. A truncated basis is a sound set of consequences of the system but need not be complete, so "not in the basis" does not mean "not a consequence".
prolongation_roundsint: Number of prolongation rounds that contributed new relations.

Properties

consistent

#source

truncated

#source

prolongation_rounds

#source

Methods

working_dae()

#source

The prolonged DAE: the input system plus the derivative jets introduced while differentiating it.

variables()

#source

The jet variables indexing the basis, in exponent-slot order.

These are the symbols the elimination actually ran over — the time variable, the declared states and derivatives, and every higher jet (d2x/dt2, …) introduced by prolongation.

final_basis()

#source

The saturated Gröbner basis, or None when the system is inconsistent.

The returned GroebnerBasis knows its variable ordering, so final_basis().to_exprs() gives the eliminated relations as Expr.

Example:

Python

r = alkahest.rosenfeld_groebner(dae, max_prolong_rounds=1)
for eq in r.final_basis().to_exprs():
    print(eq, "= 0")

ODE and DAE API

#ODE

#DAE

#x' = u, y' = w, u' = -lamx, w' = -lamy - 1, x2 + y2 = 1

#Sensitivity analysis

#Acausal modeling

#Hybrid systems

#Differential algebra

ODE

DAE

x' = u, y' = w, u' = -lamx, w' = -lamy - 1, x2 + y2 = 1

Sensitivity analysis

Acausal modeling

Hybrid systems

Differential algebra