API Reference¶

Functions¶

class Solverz.sym_algebra.functions.sin(*args)¶: The sine function.

class Solverz.sym_algebra.functions.cos(*args)¶: The cosine function.

class Solverz.sym_algebra.functions.exp(*args)¶: The exponential function, \(e^x\).

class Solverz.sym_algebra.functions.Abs(*args)¶

The element-wise absolute value function:

\[\operatorname{Abs}(x)=|x|\]

with derivative

\[\frac{\mathrm{d}}{\mathrm{d}x}{\operatorname{Abs}}(x)=\operatorname{Sign}(x).\]

class Solverz.sym_algebra.functions.Sign(*args)¶: The element-wise indication of the sign of a number

\[\begin{split}\operatorname{Sign}(x)= \begin{cases} 1&x> 0\\ 0& x==0\\ -1&x< 0 \end{cases}\end{split}\]

class Solverz.sym_algebra.functions.AntiWindUp(u, umin, umax, e)¶

For PI controller

\[\begin{split}\begin{aligned} \dot{z}(t) & =e(t) \\ u_{d e s}(t) & =K_p e(t)+K_i z(t) \end{aligned},\end{split}\]

we limit the integrator output by setting, in the \(K_i>0\) case,

\[\begin{split}\dot{z}(t)= \begin{cases} 0 & \text { if } u_{\text {des }}(t) \geq u_{\max } \text { and } e(t) \geq 0 \\ 0 & \text { if } u_{\text {des }}(t) \leq u_{\min } \text { and } e(t) \leq 0 \\ e(t) & \text { otherwise } \end{cases}.\end{split}\]

So the AntiWindUp function reads

\[\dot{z}(t)= \operatorname{AntiWindUp}(u_{\text {des }}(t), u_{\min }, u_{\max }, e(t)).\]

class Solverz.sym_algebra.functions.Min(x, y)¶: The element-wise minimum of two numbers

\[\begin{split}\begin{split}\min(x,y)= \begin{cases} x&x\leq y\\ y& x>y \end{cases}\end{split}\end{split}\]

class Solverz.sym_algebra.functions.MatVecMul(*args)¶

Legacy sparse matrix-vector multiply. Use Mat_Mul instead.

MatVecMul decomposes the sparse matrix into CSC components and calls a Numba-compatible csc_matvec. New code should use Mat_Mul(A, x) which uses scipy.sparse directly and supports full matrix calculus.

class Solverz.sym_algebra.functions.Saturation(*args)¶: The element-wise saturation a number

\[\begin{split}\operatorname{Saturation}(v, v_\min, v_\max)= \begin{cases} v_\max&v> v_\max\\ v& v_\min\leq v\leq v_\max\\ v_\min&v< v_\min \end{cases}\end{split}\]

class Solverz.sym_algebra.functions.heaviside(*args)¶

The heaviside step function

\[\begin{split}\operatorname{Heaviside}(x)= \begin{cases} 1 & x >= 0\\ 0 & x < 0\\ \end{cases}\end{split}\]

which should be distinguished from sympy.Heaviside

class Solverz.sym_algebra.functions.ln(*args)¶: The ln function, \(ln(x)\).

Equations¶

class Solverz.equation.eqn.Eqn(name: str, eqn)¶: The Equation object

class Solverz.equation.eqn.Ode(name: str, f, diff_var: iVar | IdxVar | Var)¶

The class for ODE reading

\[\frac{\mathrm{d}y}{\mathrm{d}t}=f(t,y)\]

where \(y\) is the state vector.

class Solverz.equation.eqn.LoopEqn(name: str, outer_index=None, body=None, n_outer: int = None, var_map: Dict[str, object] = None, model=None, outer=None)¶

Equation declared as a parameterised scalar template that prints to a Python for-loop in the generated inner_F / inner_J.

The motivation: SolMuseum’s network modules (eps_network, heat_network, gas_network) write per-bus / per-node equations using Python for-loops that emit 2 * n_bus (or more) scalar Eqn objects. Each scalar Eqn becomes its own inner_F<N> sub-function, and on the IES benchmark the resulting Numba LLVM compile time dominates the warm-up phase. LoopEqn expresses the same thing as a SINGLE Eqn whose body is a scalar template parameterised by an outer index i; the printer then emits ONE inner_F<N> containing a for i in range(n_outer): loop. One sub-function instead of n_outer.

Parameters:

name : str

Equation name.

outer_index : sympy.Idx

The outer loop index. Each row of the residual is the body evaluated at one value of outer_index.

n_outer : int

Range of the outer loop. The residual block has size n_outer.

body : sympy.Expr

A scalar sympy expression built either from sympy.IndexedBase (+ var_map) or from Solverz IdxVar / IdxPara (+ model). May involve the outer_index, inner sympy.Sum over auxiliary Idx symbols, and arithmetic.

model : Model, optional

When the body uses Solverz IdxVar / IdxPara, pass the containing Model so LoopEqn can resolve each name0 back to a real Var / Param (needed for sparsity detection and CSR pre-computation). Mutually exclusive with var_map.

var_map : Dict[str, sSymBasic], optional

Legacy API. Maps each sympy.IndexedBase name in body to the corresponding Solverz Var / Param instance. Mutually exclusive with model.

Substrate

Internally, the body is rewritten as a sympy expression built from sympy.IndexedBase, sympy.Idx, sympy.Sum. Solverz’s own IdxVar / IdxPara are opaque to sympy’s subs / Sum.doit() / diff, so we can’t use them verbatim — LoopEqn walks the body once at construction and rewrites any Solverz IdxVar / IdxPara into the corresponding sympy.IndexedBase form. The user is free to write either:

Native Solverz syntax (recommended): m.G[i, j] and m.ux[j] directly inside the body, and pass model=m so LoopEqn can look up each name in the Model to resolve sparsity / dim / value info.
Explicit sympy syntax (legacy): declare G_sym = sp.IndexedBase('G') etc., build the body out of those, and pass a var_map={'G': m.G, ...} dict mapping each IndexedBase name to a Solverz Var/Param.

Both styles produce the same internal representation. Native syntax is terser and avoids the bookkeeping trap of forgetting to add an entry to var_map when adding a new reference.

Examples

Native Solverz syntax (preferred):

m.ux = Var('ux', np.ones(nb))
m.uy = Var('uy', np.zeros(nb))
m.G = Param('G', csc_array(G_dense), dim=2, sparse=True)
m.B = Param('B', csc_array(B_dense), dim=2, sparse=True)
m.ix_pin = Param('ix_pin', ix_pin)

i, j = sp.Idx('i'), sp.Idx('j')
body = (m.ix_pin[i]
        - sp.Sum(m.G[i, j] * m.ux[j], (j, 0, nb - 1))
        + sp.Sum(m.B[i, j] * m.uy[j], (j, 0, nb - 1)))

m.ix_inj = LoopEqn(
    'ix_inj',
    outer_index=i, n_outer=nb,
    body=body,
    model=m,
)

Legacy var_map API (still supported):

ux_sym = sp.IndexedBase('ux')
G_sym = sp.IndexedBase('G')
body = sp.Sum(G_sym[i, j] * ux_sym[j], (j, 0, nb - 1))
m.ix_inj = LoopEqn(
    'ix_inj',
    outer_index=i, n_outer=nb,
    body=body,
    var_map={'ux': m.ux, 'G': m.G},
)

The generated inner_F for either form is a single function containing one Python for i in range(nb): loop, replacing the 2 * nb scalar sub-functions the original loop pattern would emit.

class Solverz.equation.eqn.LoopOde(name: str, outer_index=None, body=None, diff_var: iVar | IdxVar | Var = None, n_outer: int = None, model=None, var_map=None, outer=None)¶

A LoopEqn whose residual is the time derivative of a contiguous state-Var slice, i.e. an Ode with a scalar-template body.

\[\frac{\mathrm{d}y[i]}{\mathrm{d}t} = f(i,\, y,\, \text{params}) \quad\text{for}\ i = 0, \ldots, n_{\text{outer}} - 1\]

where y (diff_var) is a length-n_outer Var slice and the body f is a scalar sympy expression parameterised by the LoopEqn outer_index.

Compared to LoopEqn, LoopOde:

Overrides LHS to sp.Derivative(diff_var, t) so the rest of Solverz’s DAE assembly (eqn_dict partitioning into f_list vs g_list, mass-matrix construction, eval_lhs) sees it as an Ode and routes it to the differential-equation side of the state vector.
Participates in the LoopEqn module-printer dispatch via the existing isinstance(eqn, LoopEqn) branch, so it emits exactly ONE inner_F<N> sub-function with a Python for loop instead of one per state-cell (or per per-pipe Ode as the legacy heat-PDE pattern does). The residual vector returned by that sub-function is written to _F_[diff_var_state_slice] through the standard print_eqn_assignment_with_precompute path.

Multi-inheritance from LoopEqn and Ode (Eqn is the shared base) gives isinstance checks on both the LoopEqn-specific module-printer path and the Ode-specific DAE-split path the right answer. MRO puts LoopEqn before Ode, so LoopEqn.__init__ is what runs first; we then manually overwrite LHS and attach diff_var.

Parameters:

name : str

Equation name (same as LoopEqn).

outer_index : sympy.Idx

Scalar-template outer loop index.

body : sympy.Expr

Scalar body expression (same semantics as LoopEqn).

diff_var : iVar / IdxVar / Var

The state-Var slice being differentiated in time. Its runtime length must match n_outer.

model, var_map, n_outer : same as LoopEqn.

class Solverz.equation.eqn.IndexSet(name: str, values)¶

Named set of integer indices used as the outer iteration range of a LoopEqn / LoopOde.

Pyomo’s Set primitive inspired the API. Two construction forms:

Identity range — Set('Nodes', n) where n is a positive int. The set is the full range [0, n). The translator emits a plain for i in range(n): loop; no auxiliary Param is stored.
Explicit subset — Set('PVPQ', values) where values is a 1-D integer sequence (np.ndarray, list, tuple, or range — anything np.asarray(..., dtype=int) can eat). The values must be non-negative and duplicate-free. An internal Param named after the set is created lazily when the set is registered on a Model so the generated kernel can read for _i in range(len(values)): i = set_param[_i].

The set exposes four attributes the LoopEqn machinery consumes:

size — length of the set.
is_identity — True when values are [0, len) so the direct-outer (no Param) fast path can be used.
param — the auxiliary 1-D int Param for a non-identity set; None for an identity set. Lazily created when _register() is called by Model.__setattr__.
values — the underlying np.ndarray of indices.

Examples

Identity range (equivalent to Idx('i', n) + explicit n_outer):

m.Bus = Set('Bus', nb)
i = m.Bus.idx('i')
m.eqn = LoopEqn('eqn', outer=m.Bus, body=body, model=m)

Explicit subset (replaces the Param + Idx + indirect-access pattern):

m.PVPQ = Set('PVPQ', pv_pq_array)
i = m.PVPQ.idx('i')
m.eqn = LoopEqn('eqn', outer=m.PVPQ, body=m.Vm[i], model=m)

idx(name: str) → Idx¶

Produce a bounded sympy.Idx sized to this set.

Body authors use the returned Idx exactly as they would a bare Idx(name, n) — the body rewriter detects that the Idx came from a non-identity set and transparently wraps every Var / Param access with the set’s gather Param.

property param¶

The auxiliary 1-D int Param that stores the set’s values at runtime. Materialised lazily on first access so body rewriters can reference it without waiting for Model.add.

Identity full-range sets (values == arange(size)) still receive a Param — at the API layer is_identity only tells the translator “no gather needed when the target is the same size”, but small-identity sets used as subsets of a larger target (e.g. Ref = Set('Ref', [0]) indexing a 30-bus Var) still need the Param at runtime to emit the gather.

Solverz.equation.eqn.Idx(name: str, n: int = None)¶

Create a sympy sympy.Idx for LoopEqn body construction — a thin shortcut so users don’t need to import sympy as sp just to grab one index symbol.

Parameters:

name : str

The index label.

n : int, optional

If given, the index is bounded to [0, n - 1] (sympy stores this as .lower / .upper). Sum() and LoopEqn can then pick up the range automatically:
>>> i = Idx('i', 3)
>>> Sum(m.G[i, j] * m.ux[j], j)  # no explicit range
>>> m.eqn = LoopEqn('eqn', outer_index=i, body=body, model=m)
# n_outer auto-inferred from i.upper - i.lower + 1 = 3
Without n the result is an unbounded sp.Idx and Sum() / LoopEqn require explicit range / n_outer.

Solverz.equation.eqn.Sum(expr, dummy, n: int = None)¶

Range-aware shortcut for sympy.Sum used by LoopEqn bodies. Three call forms:

Sum(expr, j, n) → sp.Sum(expr, (j, 0, n - 1)) — the most common case. n is an int.
Sum(expr, j) where j = Idx('j', n) is a bounded sp.Idx → sp.Sum(expr, (j, j.lower, j.upper)).
Sum(expr, (j, lo, hi)) → passthrough to sp.Sum (legacy sympy tuple form).

Raises ValueError if neither n nor a bounded Idx is available.

Utilities¶

Solverz.utilities.miscellaneous.derive_incidence_matrix(indices, m, n) → coo_array¶: Convert indices to m*n incidence matrix

Solvers¶

AE solver¶

Solverz.solvers.nlaesolver.nr.nr_method(eqn: nAE, y: ndarray, opt: Opt = None)¶

The Newton-Raphson method.

Parameters:

eqn : nAE

Numerical AE object.

y : np.ndarray

The initial values of variables

opt : Opt

The solver options, including:

ite_tol: 1e-5(default)|float
The iteration error tolerance.

Returns:

sol : aesol

The aesol object.

References

[R1]

https://en.wikipedia.org/wiki/Newton%27s_method

Solverz.solvers.nlaesolver.cnr.continuous_nr(eqn: nAE, y: ndarray, opt: Opt = None)¶

The Continuous Newton-Raphson method.

Parameters:

eqn : nAE

Numerical AE object.

y : np.ndarray

The initial values of variables

opt : Opt

The solver options, including:

ite_tol: 1e-5(default)|float
The iteration error tolerance.

Returns:

sol : aesol

The aesol object.

References

[R2]

F. Milano, ‘Continuous Newton’s Method for Power Flow Analysis’, IEEE Transactions on Power Systems, vol. 24, no. 1, pp. 50–57, Feb. 2009, doi: 10.1109/TPWRS.2008.2004820.

Solverz.solvers.nlaesolver.lm.lm(eqn: nAE, y: ndarray, opt: Opt = None)¶

The Levenberg-Marquardt method by minpack, wrapped in scipy.

Warning

Note that this function uses only dense Jacobian.

Parameters:

eqn : nAE

Numerical AE object.

y : np.ndarray

The initial values of variables

opt : Opt

The solver options, including:

ite_tol: 1e-5(default)|float
The iteration error tolerance.

Returns:

sol : aesol

The aesol object.

References

[R3]

More, Jorge J., Burton S. Garbow, and Kenneth E. Hillstrom. 1980. User Guide for MINPACK-1.

[R4]

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.root.html

Solverz.solvers.nlaesolver.sicnm.sicnm(ae: nAE, y0: ndarray, opt: Opt = None)¶

The semi-implicit continuous Newton method [R5]. The philosophy is to rewrite the algebraic equation

\[0=g(y)\]

as the differential algebraic equations

\[\begin{split}\begin{aligned} \dot{y}&=z \\ 0&=J(y)z+g(y) \end{aligned}\end{split}\]

with \(y_0\) being the initial value guess and \(z_0=-J(y_0)^{-1}g(y_0)\), where \(z\) is an intermediate variable introduced. Then the DAEs are solved by Rodas. SICNM is found to be more robust than the Newton’s method, for which the theoretical proof can be found in my paper [R5]. In addition, the non-iterative nature of Rodas guarantees the efficiency.

One can change the rodas scheme according to the ones implemented in the DAE version of Rodas.

SICNM used the Hessian-vector product (HVP) interface of the algebraic equation models, the make_hvp flag should be set to True when printing numerical modules.

An illustrative example of SICNM can be found in the power flow section of Solverz’ cookbook.

Parameters:

ae : nAE

Numerical AE object.

y0 : np.ndarray

The initial values of variables

opt : Opt

The solver options, including:

ite_tol: 1e-5(default)|float
The iteration error tolerance.

Returns:

sol : aesol

The aesol object.

References

[R5] (1,2,3)

Yu, W. Gu, S. Lu, and Y. Xu, “Semi-implicit continuous newton method for power flow analysis,” 2023, arXiv:2312.02809.

FDAE solver¶

Solverz.solvers.fdesolver.fdae_solver(fdae: nFDAE, tspan: List | ndarray, y0: ndarray, opt: Opt = None, **kwargs)¶

The general solver of FDAE.

Parameters:

fdae : nFDAE

Numerical FDAE object.

tspan : List | np.ndarray

An array specifying t0 and tend

y0 : np.ndarray

The initial values of variables

opt : Opt

The solver options, including:

step_size: 1e-3(default)|float
The step size

ite_tol: 1e-8(default)|float
The error tolerance of inner Newton iterations.

kwargs : dict

The value of y at previous time nodes. For example, if an FDAE uses \(y_0\), \(y_{-1}\) and \(y_{-2}\), then we can call FDAE solver with
sol = fdae_solver(mdl,
                  [0, 120],
                  y0,
                  Opt(step_size=60),
                  y1=y1,
                  y2=y2)
where y1 denotes \(y_{-1}\) and y2 denotes \(y_{-2}\).

Returns:

sol : daesol

The daesol object.

DAE solver¶

Solverz.solvers.daesolver.beuler.backward_euler(dae: nDAE, tspan: List | ndarray, y0: ndarray, opt: Opt = None)¶

The fixed step backward Euler (implicit Euler) method.

Parameters:

dae : nDAE

Numerical DAE object.

tspan : List | np.ndarray

An array specifying t0 and tend

y0 : np.ndarray

The initial values of variables

opt : Opt

The solver options, including:

step_size: 1e-3(default)|float
The step size

Returns:

sol : daesol

The daesol object.

References

[R6]

https://en.wikipedia.org/wiki/Backward_Euler_method

Solverz.solvers.daesolver.trapezoidal.implicit_trapezoid(dae: nDAE, tspan: List | ndarray, y0: ndarray, opt: Opt = None)¶

The fixed step implicit trapezoidal method.

Parameters:

dae : nDAE

Numerical DAE object.

tspan : List | np.ndarray

An array specifying t0 and tend

y0 : np.ndarray

The initial values of variables

opt : Opt

The solver options, including:

step_size: 1e-3(default)|float
The step size

Returns:

sol : daesol

The daesol object.

References

[R7]

https://en.wikipedia.org/wiki/Trapezoidal_rule_(differential_equations)

Solverz.solvers.daesolver.rodas.Rodas(dae: nDAE, tspan: List | ndarray, y0: ndarray, opt: Opt = None)¶

The stiffly accurate Rosenbrock methods including Rodas4 [R8], Rodasp [R9], Rodas5p [R10].

Parameters:

dae : nDAE

Numerical DAE object.

tspan : List | np.ndarray

Either - a list specifying [t0, tend], or - a np.ndarray specifying the time nodes that you are concerned about

y0 : np.ndarray

The initial values of variables

opt : Opt

The solver options, including:

scheme: ‘rodas4’(default)|’rodasp’|’rodas5p’
The rodas scheme

rtol: 1e-3(default)|float
The relative error tolerance

atol: 1e-6(default)|float
The absolute error tolerance

event: Callable
The simulation events, with \(t\) and \(y\) being args, and value, is_terminal and direction being outputs

fixed_h: False(default)|bool
To use fixed step size

hinit: None(default)|float
Initial step size

hmax: None(default)|float
Maximum step size.

pbar: False(default)|bool
To display progress bar

Returns:

sol : daesol

The daesol object.

References

[R8] (1,2)

Hairer and Wanner, Solving Ordinary Differential Equations II , 2nd ed. Berlin, Heidelberg, Germany: Springer-Verlag, 1996.

[R9] (1,2)

Steinebach, Order-reduction of ROW-methods for DAEs and method of lines applications. Preprint-Nr. 1741, FB Mathematik, TH Darmstadt, 1995

[R10] (1,2)

Steinebach, “Construction of rosenbrock-wanner method rodas5p and numerical benchmarks within the julia differential equations package,” BIT, vol. 63, no. 27, Jun 2023.

Solverz.solvers.daesolver.radau.Radau(dae: nDAE, tspan: List | ndarray, y0: ndarray, opt: Opt = None)¶

Radau IIA(5) DAE solver — 3-stage 5th-order fully-implicit Runge-Kutta method, ported from SciML’s RadauIIA5.

Parameters:

dae : nDAE

Numerical DAE object.

tspan : List | np.ndarray

Either - a list specifying [t0, tend], or - a np.ndarray specifying the time nodes that you are concerned about

y0 : np.ndarray

The initial values of variables

opt : Opt

The solver options, including:

rtol: 1e-3(default)|float
The relative error tolerance

atol: 1e-6(default)|float
The absolute error tolerance

event: Callable
The simulation events, with \(t\) and \(y\) being args, and value, is_terminal and direction being outputs

hinit: None(default)|float
Initial step size

hmax: None(default)|float
Maximum step size

max_it: 100(default)|int
The simplified Newton maximum iteration count per stage attempt. SciML’s Radau5 default is 10; the lower bound is applied automatically.

pbar: False(default)|bool
To display progress bar

Returns:

sol : daesol

The daesol object.

Notes

Stiffly accurate: y_{n+1} = y_n + Z_3. The Newton update is performed in the W-coordinate system (W = TI @ Z), which decouples the linear solves into one real and one complex sparse LU factorisation per step attempt. Step-size adaptation uses the SciML/Hairer predictive (Gustafsson) controller.

References

[R11]

E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed. Berlin, Heidelberg: Springer-Verlag, 1996, §IV.8.

[R12]

SciML/OrdinaryDiffEq.jl, OrdinaryDiffEqFIRK module, RadauIIA5.

Solverz.solvers.daesolver.ode15s.ode15s(dae: nDAE, tspan: List | ndarray, y0: ndarray, opt: Opt = None)¶

The python implementation of MATLAB Ode15s without event detection [R13].

Parameters:

dae : nDAE

Numerical DAE object.

tspan : List | np.ndarray

Either - a list specifying [t0, tend], or - a np.ndarray specifying the time nodes that you are concerned about

y0 : np.ndarray

The initial values of variables

opt : Opt

The solver options, including:

rtol: 1e-3(default)|float
The relative error tolerance

atol: 1e-6(default)|float
The absolute error tolerance

hinit: None(default)|float
Initial step size

hmax: None(default)|float
Maximum step size.

pbar: False(default)|bool
To display progress bar

Returns:

sol : daesol

The daesol object.

References

[R13] (1,2)

1. Shampine and M. W. Reichelt, “The MATLAB ODE Suite,” SIAM J. Sci. Comput., vol. 18, no. 1, pp. 1–22, Jan. 1997, doi: 10.1137/S1064827594276424.