Implementation of Numerical integration, equilibrium finding, and discretization utilities. More...

#include "Differential_equations.h"
#include <stdexcept>
#include <iostream>
#include <cmath>

Include dependency graph for Differential_equations.cpp:

Classes
struct	CVODEUserData

struct	KINSOLUserData

Functions
static Eigen::VectorXd	nv2eigen (N_Vector v, int n)

static void	eigen2nv (N_Vector v, const Eigen::VectorXd &x)

static void	eigen2sun (SUNMatrix J, const Eigen::MatrixXd &M)

static int	cvode_rhs_cb (sunrealtype t, N_Vector y, N_Vector ydot, void *ud)

static int	cvode_jac_cb (sunrealtype t, N_Vector y, N_Vector, SUNMatrix J, void *ud, N_Vector, N_Vector, N_Vector)

static int	kinsol_res_cb (N_Vector u, N_Vector fval, void *ud)

static int	kinsol_jac_cb (N_Vector u, N_Vector, SUNMatrix J, void *ud, N_Vector, N_Vector)

static int	kinsol_fd_jac_cb (N_Vector u, N_Vector, SUNMatrix J, void *ud, N_Vector, N_Vector)

static KINSOLConfig	withAutoScaling (const KINSOLConfig &cfg, const RHSFunc &rhs, const Eigen::VectorXd &x0, const Eigen::VectorXd &u)

static Eigen::VectorXd	kinsolSolve (const RHSFunc &rhs, const Eigen::VectorXd &x0, const Eigen::VectorXd &u, const KINSOLConfig &cfg, const JacFunc &jac)

Eigen::VectorXd	findEquilibrium (const RHSFunc &rhs, const Eigen::VectorXd &x0, const Eigen::VectorXd &u, const KINSOLConfig &cfg, const JacFunc &jac)
	Finds a steady-state equilibrium dx/dt = 0 for fixed input u.

Eigen::VectorXd	findEquilibriumRobust (const RHSFunc &rhs, const Eigen::VectorXd &x0, const Eigen::VectorXd &u, const JacFunc &jac)
	Robust equilibrium finder with a cascade of solver strategies.

CVODEResult	integrate (const RHSFunc &rhs, const Eigen::VectorXd &x0, const std::function< Eigen::VectorXd(double t)> &inputFn, double t0, double tEnd, double dt_output, const CVODEConfig &cfg, const JacFunc &jac)
	Integrates an ODE system using CVODE (BDF, stiff).

std::pair< Eigen::MatrixXd, Eigen::MatrixXd >	computeJacobians (const RHSFunc &rhs, const Eigen::VectorXd &x, const Eigen::VectorXd &u, double t, double eps)
	Computes ∂f/∂x and ∂f/∂u by central finite differences.

void	padeDelaySystem3 (double T, MatrixXd &Ad, MatrixXd &Bd, MatrixXd &Cd, MatrixXd &Dd)

void	padeDelaySystemMulti3 (double T, MatrixXd &Ad, MatrixXd &Bd, MatrixXd &Cd, MatrixXd &Dd, int ns)

void	padeDelaySystem2 (double T, MatrixXd &Ad, MatrixXd &Bd, MatrixXd &Cd, MatrixXd &Dd)

void	padeDelaySystemMulti2 (double T, MatrixXd &Ad, MatrixXd &Bd, MatrixXd &Cd, MatrixXd &Dd, int ns)

void	discretizeABCD (const Eigen::MatrixXd &A, const Eigen::MatrixXd &B, const Eigen::MatrixXd &C, const Eigen::MatrixXd &D, double Ts, Eigen::MatrixXd &Ad, Eigen::MatrixXd &Bd, Eigen::MatrixXd &Cd, Eigen::MatrixXd &Dd)
	Discretizes a continuous-time state-space model using the Tustin (bilinear) method.

Detailed Description

Implementation of Numerical integration, equilibrium finding, and discretization utilities.

Function Documentation

◆ computeJacobians()

std::pair< Eigen::MatrixXd, Eigen::MatrixXd > computeJacobians	(	const RHSFunc &	rhs,
		const Eigen::VectorXd &	x,
		const Eigen::VectorXd &	u,
		double	t = `0.0`,
		double	eps = `1e-8`
	)

Computes ∂f/∂x and ∂f/∂u by central finite differences.

Parameters

rhs	ODE right-hand side.
x	State at which to evaluate.
u	Input at which to evaluate.
t	Evaluation time (default 0).
eps	Finite-difference step size (default 1e-8).

Returns: Pair (df/dx, df/du).

References j.

Referenced by MMC::computeABCD(), MMC::computeABCD_analytical(), and kinsol_fd_jac_cb().

◆ cvode_jac_cb()

static int cvode_jac_cb	(	sunrealtype	t,
		N_Vector	y,
		N_Vector	,
		SUNMatrix	J,
		void *	ud,
		N_Vector	,
		N_Vector	,
		N_Vector
	)

static

References eigen2sun(), and nv2eigen().

Referenced by integrate().

◆ cvode_rhs_cb()

static int cvode_rhs_cb	(	sunrealtype	t,
		N_Vector	y,
		N_Vector	ydot,
		void *	ud
	)

static

References eigen2nv(), and nv2eigen().

Referenced by integrate().

◆ discretizeABCD()

void discretizeABCD	(	const Eigen::MatrixXd &	A,
		const Eigen::MatrixXd &	B,
		const Eigen::MatrixXd &	C,
		const Eigen::MatrixXd &	D,
		double	Ts,
		Eigen::MatrixXd &	Ad,
		Eigen::MatrixXd &	Bd,
		Eigen::MatrixXd &	Cd,
		Eigen::MatrixXd &	Dd
	)

Discretizes a continuous-time state-space model using the Tustin (bilinear) method.

Parameters

A	Continuous-time state matrix.
B	Continuous-time input matrix.
C	Continuous-time output matrix.
D	Continuous-time feed-through matrix.
Ts	Sample period (s).
Ad	Output discrete-time state matrix.
Bd	Output discrete-time input matrix.
Cd	Output discrete-time output matrix.
Dd	Output discrete-time feed-through matrix.

Referenced by Converter::discretize(), and DQsym::run().

◆ eigen2nv()

static void eigen2nv	(	N_Vector	v,
		const Eigen::VectorXd &	x
	)

static

Referenced by cvode_rhs_cb(), integrate(), kinsol_res_cb(), and kinsolSolve().

◆ eigen2sun()

static void eigen2sun	(	SUNMatrix	J,
		const Eigen::MatrixXd &	M
	)

static

References j.

Referenced by cvode_jac_cb(), kinsol_fd_jac_cb(), and kinsol_jac_cb().

◆ findEquilibrium()

Eigen::VectorXd findEquilibrium	(	const RHSFunc &	rhs,
		const Eigen::VectorXd &	x0,
		const Eigen::VectorXd &	u,
		const KINSOLConfig &	cfg = `{}`,
		const JacFunc &	jac = `nullptr`
	)

Finds a steady-state equilibrium dx/dt = 0 for fixed input u.

Parameters

rhs	ODE right-hand side.
x0	Initial guess for the state vector.
u	Fixed input vector.
cfg	Solver strategy and tolerance settings.
jac	Optional analytical Jacobian.

Returns: Equilibrium state vector.

References kinsolSolve().

◆ findEquilibriumRobust()

Eigen::VectorXd findEquilibriumRobust	(	const RHSFunc &	rhs,
		const Eigen::VectorXd &	x0,
		const Eigen::VectorXd &	u,
		const JacFunc &	jac = `nullptr`
	)

Robust equilibrium finder with a cascade of solver strategies.

Attempts LineSearch, then Newton, relaxed warmup, and Picard + Newton until convergence or all strategies are exhausted.

Parameters

rhs	ODE right-hand side.
x0	Initial guess for the state vector.
u	Fixed input vector.
jac	Optional analytical Jacobian.

Returns: Equilibrium state vector.

References KINSOLConfig::damping, KINSOLConfig::ftol, kinsolSolve(), LineSearch, KINSOLConfig::max_iter, Newton, Picard, KINSOLConfig::stol, KINSOLConfig::strategy, and withAutoScaling().

Referenced by MMC::solveEquilibrium().

◆ integrate()

CVODEResult integrate	(	const RHSFunc &	rhs,
		const Eigen::VectorXd &	x0,
		const std::function< Eigen::VectorXd(double t)> &	inputFn,
		double	t0,
		double	tEnd,
		double	dt_output,
		const CVODEConfig &	cfg = `{}`,
		const JacFunc &	jac = `nullptr`
	)

Integrates an ODE system using CVODE (BDF, stiff).

Parameters

rhs	Right-hand side f(t, x, u).
x0	Initial state vector.
inputFn	Time-varying input u(t).
t0	Start time.
tEnd	End time.
dt_output	Output sampling interval.
cfg	Tolerance and Jacobian settings.
jac	Optional analytical Jacobian (nullptr → finite differences).

Returns: Sampled time and state histories.

References CVODEConfig::atol, cvode_jac_cb(), cvode_rhs_cb(), CVODEConfig::dt_max, eigen2nv(), nv2eigen(), CVODEConfig::rtol, CVODEResult::states, CVODEResult::time, and CVODEConfig::use_analytical_jac.

◆ kinsol_fd_jac_cb()

static int kinsol_fd_jac_cb	(	N_Vector	u,
		N_Vector	,
		SUNMatrix	J,
		void *	ud,
		N_Vector	,
		N_Vector
	)

static

References computeJacobians(), eigen2sun(), and nv2eigen().

Referenced by kinsolSolve().

◆ kinsol_jac_cb()

static int kinsol_jac_cb	(	N_Vector	u,
		N_Vector	,
		SUNMatrix	J,
		void *	ud,
		N_Vector	,
		N_Vector
	)

static

References eigen2sun(), and nv2eigen().

Referenced by kinsolSolve().

◆ kinsol_res_cb()

static int kinsol_res_cb	(	N_Vector	u,
		N_Vector	fval,
		void *	ud
	)

static

References eigen2nv(), and nv2eigen().

Referenced by kinsolSolve().

◆ kinsolSolve()

static Eigen::VectorXd kinsolSolve	(	const RHSFunc &	rhs,
		const Eigen::VectorXd &	x0,
		const Eigen::VectorXd &	u,
		const KINSOLConfig &	cfg,
		const JacFunc &	jac
	)

static

◆ nv2eigen()

static Eigen::VectorXd nv2eigen	(	N_Vector	v,
		int	n
	)

static

Referenced by cvode_jac_cb(), cvode_rhs_cb(), integrate(), kinsol_fd_jac_cb(), kinsol_jac_cb(), kinsol_res_cb(), and kinsolSolve().

◆ padeDelaySystem2()

void padeDelaySystem2	(	double	T,
		MatrixXd &	Ad,
		MatrixXd &	Bd,
		MatrixXd &	Cd,
		MatrixXd &	Dd
	)

◆ padeDelaySystem3()

void padeDelaySystem3	(	double	T,
		MatrixXd &	Ad,
		MatrixXd &	Bd,
		MatrixXd &	Cd,
		MatrixXd &	Dd
	)

◆ padeDelaySystemMulti2()

void padeDelaySystemMulti2	(	double	T,
		MatrixXd &	Ad,
		MatrixXd &	Bd,
		MatrixXd &	Cd,
		MatrixXd &	Dd,
		int	ns
	)

Referenced by MMC::MMC(), and MMC::MMC().

◆ padeDelaySystemMulti3()

void padeDelaySystemMulti3	(	double	T,
		MatrixXd &	Ad,
		MatrixXd &	Bd,
		MatrixXd &	Cd,
		MatrixXd &	Dd,
		int	ns
	)

Referenced by MMC::MMC(), and MMC::MMC().

◆ withAutoScaling()

static KINSOLConfig withAutoScaling	(	const KINSOLConfig &	cfg,
		const RHSFunc &	rhs,
		const Eigen::VectorXd &	x0,
		const Eigen::VectorXd &	u
	)

static

References KINSOLConfig::f_scale, and KINSOLConfig::x_scale.

Referenced by findEquilibriumRobust().

Classes

Functions

Detailed Description

Function Documentation

◆ computeJacobians()

◆ cvode_jac_cb()

◆ cvode_rhs_cb()

◆ discretizeABCD()

◆ eigen2nv()

◆ eigen2sun()

◆ findEquilibrium()

◆ findEquilibriumRobust()

◆ integrate()

◆ kinsol_fd_jac_cb()

◆ kinsol_jac_cb()

◆ kinsol_res_cb()

◆ kinsolSolve()

◆ nv2eigen()

◆ padeDelaySystem2()

◆ padeDelaySystem3()

◆ padeDelaySystemMulti2()

◆ padeDelaySystemMulti3()

◆ withAutoScaling()