Getting Started¶

Overview¶

Solverz supports three types of equality constraints, which are respectively

Algebraic Equations (AEs) with general formula \(0=F(y,p)\)
Finite Difference Algebraic Equations (FDAEs) with general formula \(0=F(t,y,p,y_0)\)
Differential Algebraic Equations (DAEs) with general formula \(M\dot{y}=F(t,y,p)\)

where \(y\) is the known variable, \(F\) is the mapping, \(p\) is the parameter set of your models, \(y_0\) is the previous time node value of \(y\), \(M\)is the singular mass matrix.

If your mathematical models belong to one of the equation types above, then Solverz should be your choice, because it supports an object-oriented design for the definition of these equations and can generate high-performance numerical codes.

The basic steps of a simple modeling and solution process are:

Create model and declare components
Instantiate the model and generate numerical codes
Apply solver
Interrogate solver results

These steps and the philosophy behind Solverz are explained as follows.

Symbolic Modelling¶

The modelling starts with declaring an empty Model with

from Solverz import Model
m = Model()

Variables¶

Variables represent unknown or changing parts of a model. The values taken by the variables are often referred to as a solution of the simulation. In Solverz, variables can be added to the model by

from Solverz import Var
m.x = Var(name='x', value=[0, 1])
m.y = Var(name='y', value=0)

We should assign two attributes to each variable, which are name and value respectively. The name is how variable m.x is represented in expressions and the value is the initial conditions of the m.x.

Solverz supports only one-dimensional variables and the length of the variables is implicitly dependent on the length of the value you give.

Index¶

Sometimes we want to access certain elements in a variable, for example, the first element of x, and the first to third elements of y. Solverz provides you with the flexibility by just calling

>>> m.x[0]
x[0]
>>> m.y[0:3]
y[0:3]

Currently, only int and slice indices are recommended. We are working on more complex index types such as the list. Also, be careful not to exceed the maximum index that depends on the initial value you give. Otherwise, errors would be raised.

Operation¶

The symbolic mathematics of Solverz variables are based on Sympy. That is, we can use the Solverz variables to perform the python operations +, -, *, /, **, just to name a few. For example, the expressions \(x_0x_1+\sin(x_2)\) can be implemented by

>>> from Solverz import Var, sin
>>> m.x = Var(name='x', value=[0, 0, 0])
>>> m.x[0]*m.x[1]+sin(m.x[2])
x[0]*x[1] + sin(x[2])

Currently, functions like sin, cos, Abs, Sign, exp, just to name a few, are supported. You can refer to api reference for all supported functions and their detailed implementations. If you want to write your own functions, refer to advanced usage. And you can contact us so that we can make it a Solverz built-in funtion.

A special case of Variables is the AliasVar, that is, the alias variables. In Solverz, the AliasVar is used to denote the historical value of some variable, which is useful in finite difference equations. For example, one can use the following codes to declare AliasVar.

from Solverz import AliasVar
m.p = Var(name='p', value=0)
m.p0 = AliasVar(name='p', init=m.p)

The name of AliasVar is used to denote that m.p0 is the historical value of m.p. And it is initialized by m.p. One can initialize the variables by assigning expressions to the init attribute, apart from giving the initial values.

Parameters¶

Parameters represents the data that must be supplied to perform the simulation. Though we can use python expression c*m.x, where c is a number, to model the expression \(cx\). But it is recommended to declare c as a new parameter if you want to change \(c\) in simulations for different results. The declaration of parameters below is the same as variable declaration.

from Solverz import Param
m.c = Param(name='c', value=0)

The parameters can also be incorporated in operations. And it can be indexed just as the variables. In Solverz, it is easy to alter the parameters in simulation, which will be shown later.

Parameters changing with time¶

Some parameters are time-varying. For example, when simulating differential equations and the finite difference algebraic equations, you want to know the results under given boundaries. In Solverz, this is realized by using TimeSeriesParam. For example, we want the parameter pb to increase from 10 to 11 in the time interval [1,2]. We can do the following.

from Solverz import TimeSeriesParam
m.pb = TimeSeriesParam('pb',
                       v_series=[10, 10, 11, 11],
                       time_series=[0, 1, 2, 100])

The v_series and time_series specify the values and time nodes of the Time-series parameter. Solverz performs linear interpolation for the value of pb at any given \(t\). If \(t>100\), pb is assumed to be 11.

Matrix parameter¶

The vertex-edge (or node-pipe) connection relationships in energy networks can be typically described by incidence matrices, which are sparse in a lot of occasions. Since version 0.6.0, Solverz begins to support a new modeling element, the two-dim matrix parameter, which can help save the modeling overhead.

An illustrative example for equation A@x-b=0, where we denote by @ the matrix-vector multiplication, is

from Solverz import Var, Param, Mat_Mul, Model, Eqn
m = Model()
m.A = Param('A', [[1, 0], [0, 1]], sparse=True, dim=2)
m.x = Var('x', [1, 2])
m.b = Param('b', [0, 0])
m.f = Eqn('f', Mat_Mul(m.A, m.x)-m.b)

Using Mat_Mul, we can write matrix-vector equations naturally. The matrix parameter A should be declared with sparse=True and dim=2. Solverz automatically computes symbolic Jacobians for these equations using matrix calculus.

Note

Solverz supports general mixed matrix-vector expressions with automatic symbolic differentiation. Operations like exp(A@x), sin(A@x), Diag(x)@A@y, transpose(A)@x, and element-wise combinations are all supported. See Matrix-Vector Calculus for details and examples.

The sparse matrix parameter A is not mutable after declaration. Its sparsity pattern is used to construct efficient Jacobian blocks.

Equations¶

Equations describe the constraints of variables. In Solverz, equations are classified into two categories: The basic equations \(0=f(y,p),\) and the ordinary differential equation (ODE) \(\dot{y}=f(y,p)\) where \(p\) denotes the parameters and \(y\) denotes the variables.

The basic equations are imposed by the Eqn objects. Below is the example of a Eqn declaration.

from Solverz import Var, Param, Eqn, made_numerical, sin, Model, nr_method

m = Model()
m.x = Var('x', 0)
m.t = Param('t', 0.3)
m.x1 = Param('x1', -0.1)
m.f = Eqn(name='f', eqn=m.x - sin(3.14 * m.x) * m.t - m.x1)

The eqn attribute should be the right hand size of \(0=f(y,p)\).

The ODEs are imposed by the Ode object. Below is the example of an ODE declaration.

from Solverz import Model, Var, Ode

m = Model()
m.h = Var('h', 0)
m.v = Var('v', 20)
m.f1 = Ode('f1', f=m.v, diff_var=m.h)

The f attribute of Ode objects denote the right hand side of \(\dot{y}=f(y,p)\) and the diff_var attribute denotes the left hand side variable.

Create Model Instance¶

After adding all the elements to the Model object, we can create the symbolic model and the numerical variable instances by calling

eqs, y0 = m.create_instance()

If you want to control the assembled order of equations or variables, you can also write

eqs, y0 = m.create_instance(eqn_sequence=['g', 'f1', 'f2'],
                            var_sequence=['v', 'h'])

This is useful when you want a particular block layout in the Jacobian.

Symbolic Equations¶

The eqs object is the symbolic equation, which can be either AE, FDAE or DAE. The type depends on the equations and variables you declare. The Solverz detects the equation and variable types, and then constructs AE, FDAE or DAE automatically. The object stores all the symbolic equations, the derivatives and the equation addresses. Also, it can be used to check if the equation number equals the variable number, which is critical for a model to be solved.

Numerical Variables¶

The y0 object is the instance of numerical variables. In Solverz, it holds the meaning of the combination of all symbolic variables Var. And it is used directly in simulation solutions. For example, if y0 contains variable h and v. Then we can access these variables by calling y0['h'] and y0['v']. If you want to set different values of h and v, you can use y0['h']=np.array([1,2,3]). It should be noted that the new values must have the same shape compared with the original ones. This is because Solverz has already allocated the addresses for the variables when creating instances. Moreover, by overloading operators, y0 can be used for addition, subtraction, multiplication, and division.

From Symbolic to Numerical¶

Numerical Equations¶

To directly use the symbolic equations for computation is too slow. Alternatively, you should use the numerical equations derived by Solverz, which are optimized for efficient simulation. Specifically, Solverz prints all the symbolic expressions to well-organized python functions based on mature libraries such as numpy, scipy and numba.

To derive numerical functions, just use the made_numerical function as follows. We have equation \(0=x-t\sin(\pi x)-x_1,\) where \(x\) is the known variabe, \(t\) and \(x_1\) are parameters. With Solverz, one has

import numpy as np
from Solverz import Var, Param, Eqn, made_numerical, sin, Model, nr_method, module_printer

m = Model()
m.x = Var('x', 0)
m.t = Param('t', 0.3)
m.x1 = Param('x1', -0.1)
m.f = Eqn('f', m.x - sin(np.pi * m.x) * m.t - m.x1)
sae, y0 = m.create_instance()
ae, code = made_numerical(sae, y0, output_code=True, sparse=True)

The ae object is a instance of the Solverz.nAE class. It omit all the symbolic expressions and other redundant details, which has three attributes:

ae.F that returns \(F(y,p)\) if one calls ae.F(y,p)
ae.J that returns the jacobian\(J(y,p)\) if one calls ae.J(y,p)
ae.p that is a python dict object storing all the parameter values.

In the above example, one has

>>> ae.F(y0, ae.p)
array([0.1])
>>> ae.J(y0, ae.p).toarray()
array([[0.0575222]])

One can alter the parameter values, for example to set \(t=1\), by

>>> ae.p['t']=np.array([1.0])

You can also inspect the numerical codes generated by Solverz with

>>> print(code['F'])
def F_(y_, p_):
    x = y_[0:1]
    t = p_["t"]
    x1 = p_["x1"]
    _F_ = zeros((1, ))
    _F_[0:1] = -t*sin(3.14159265358979*x) - x1 + x
    return _F_
>>> print(code['J'])
def J_(y_, p_):
    x = y_[0:1]
    t = p_["t"]
    x1 = p_["x1"]
    row = []
    col = []
    data = []
    row.extend([0])
    col.extend(arange(0, 1))
    data.extend(((-3.14159265358979*t*cos(3.14159265358979*x) + 1).tolist()))
    return coo_array((data, (row,col)), (1, 1)).tocsc()

Similarly, Solverz.nFDAE and Solverz.nDAE are respectively the numerical equation abstraction of FDAE and DAE, with the F and J attributes being the numerical interfaces. The Solverz.nFDAE instance has a nstep attribure to denote the number of historical time steps that is required. Currently, nstep can only be one.

If you only need a contiguous block of the symbolic Jacobian, you can call

sub_jac = eqs.FormPartialJac(y0,
                             eqn_list=['g'],
                             var_list=['v'])

The selected equation names and variable names must each be contiguous in the current ordering. The returned Jac stores the block size in shape and the global offset of the upper-left corner in coordinate0.

Sometimes, one wants to use the second derivative information. Since Release/0.1, Solverz is able to derive the Hessian-vector product, with formula

\[\frac{\partial }{\partial y}\left(J(y)z\right)=H(y)\otimes z\]

where the Hessian tensor

\[H(y)=\left(\frac{\partial \nabla g_i(y)^\mathrm{T}}{y_j}\right)_{ij}\]

and \(J(y)\) is the original Jacobian.

One can call the HVP() method of numerical equations after generating numerical modules with make_hvp=True. A successful attempt of using HVP to improve the robustness of AE’s solution can be found in Solverz’ cookbook.

Sparse matrix¶

In most cases, the Jacobian is a sparse matrix in which most of the elements are zero. The sparse matrices can be decomposed very efficiently using a sparse solver. It is suggested not to derive dense Jacobians because the decomposition and storage of zero elements in RAM can be very time-consuming.

However, Solverz still provides both sparse and dense Jacobian options. One can speficy whether to derive sparse Jacobian by setting the sparse argument in made_numerical().

The Module Printer¶

To model very complex systems can be slow. If one wants to avoid the repeated derivations of symbolic and numerical equations, a good way is to render an independent python module for your simulation model using module_printer. An illustrative example is

from Solverz import Model, Var, Ode, Opt, made_numerical, Rodas, module_printer

m = Model()
m.x = Var('x', [0, 20])
m.f1 = Ode('f1', m.x[1], m.x[0])
m.f2 = Ode('f2', -9.8, m.x[1])
bball, y0 = m.create_instance()

pyprinter = module_printer(mdl=bball,
                           variables=y0,
                           name='bounceball',
                           jit=True)
pyprinter.render()

One can import the module to the codes by

from bounceball import mdl as nbball, y as y0

The modules rendered are better optimized compared with the ones generated by made_numerical. In efficiency-demanding scenarios, we recommend you using the module printer.

Dynamic Compilation¶

The dynamic nature of Python, an interpreted-type language, brings about the efficiency decreasing compared with languages like C/C++, especially in numerical computation. To resolve this, Solverz uses Numba to accelerate your codes based on just-in-time (JIT) compilation. Specifically, Numba translates Python functions to optimized machine code at runtime using the industry-standard LLVM compiler library. Numba-compiled numerical algorithms in Python can approach the speeds of C or FORTRAN.

If one wants to take advantage of the dynamic compilation, just set the jit arg in module_printer to be True. Then the models will be compiled at the first time of function evaluation. Though the compilation time of complex models can be of tens of minutes, the compilation results are cached locally. Hence, the model should be compiled only once. It is recommended to first set jit=False to debug the models and then perform dynamic compilation. Currently, dynamic compilation is not supported in made_numerical.

Solvers¶

Solverz provides basic solvers for the solutions of AE, FDAE and DAE. We are working hard on implementing more mature solvers. Please feel free to contact us if you have any good idea~

Below is an overview of the built-in solvers.

AE solvers¶

nr_method() the Newton-Raphson method
continuous_nr() the continuous Newton method, which is more robust compared with the Newton-Raphson
lm() the Levenberg-Marquardt method provided by scipy.optimize.root. Only dense Jacobian is allowed, which may be time-consuming.
sicnm() the semi-implicit version of continuous Newton method. It possesses both the implicit stability and explicit computation overhead, which shows both robustness and efficiency. Please make Hessian-vector product if you want to use it.

FDAE solver¶

fdae_solver()

DAE solvers¶

backward_euler() the backward Euler method.
implicit_trapezoid() the implicit trapezoidal method.
Rodas() the stiffly accurate Rosenbrock method with adaptive step size, dense output and event detection. One can use it the same as the Ode-series solvers in Matlab. This is the most stable solver in Solverz.
Radau() the 3-stage 5th-order Radau IIA fully-implicit Runge-Kutta method, ported from SciML’s RadauIIA5. Stiffly accurate; recommended for stiff or oscillatory DAEs that need higher order than Rodas.

The detailed usage of these solvers can be found in api reference.

It also a good idea to use solvers provided by scipy and other python packages since Solverz has derived the generic numerical interfaces.

All solver results carry a stats object, for example sol.stats.nstep, sol.stats.nfeval, sol.stats.nJeval, sol.stats.ndecomp, sol.stats.nsolve and sol.stats.nreject. For DAE solvers, setting Opt(profile=True) prints the elapsed wall-clock time of the solver call.