Module: solver
The solver module provides interfaces to static, dynamic and eigenvalue solvers. Most of the solvers are implemented inside the C++ core.
Author: Johannes Gerstmayr
Date: 2020-12-02
Notes: Solver functions are included directly in exudyn and can be used with mbs.SolveStatic(…)
Function: SolverErrorMessage
function description: (internal) helper function for unique error and helper messages
Function: SolveStatic
function description: solves the static mbs problem using simulationSettings; see MainSolverStatic, MainSolverStatic, for further details of the static solver; this function is also available in exudyn (using exudyn.SolveStatic(…))
input:
mbs: the MainSystem containing the assembled system; note that mbs may be changed upon several runs of this functionsimulationSettings: specific simulation settings out of exu.SimulationSettings(), as described in SolutionSettings; use options for newton, discontinuous settings, etc., from staticSolver sub-itemsupdateInitialValues: if True, the results are written to initial values, such at a consecutive simulation uses the results of this simulation as the initial values of the next simulationstoreSolver: if True, the staticSolver object is stored in the mbs.sys dictionary as mbs.sys[‘staticSolver’], and simulationSettings are stored as mbs.sys[‘simulationSettings’]showHints: show additional hints, if solver failsshowCausingItems: if linear solver fails, this option helps to identify objects, etc. which are related to a singularity in the linearized system matrixautoAssemble: if True: if mbs.systemIsConsistent=False (system is not assembled), call mbs.Assemble() before solver calls
output: (type: bool) returns True, if successful, False if fails; if storeSolver = True, mbs.sys contains staticSolver, which allows to investigate solver problems (see Solver substructures and the items described in MainSolverStatic)
example:
import exudyn as exu
from exudyn.itemInterface import *
SC = exu.SystemContainer()
mbs = SC.AddSystem()
#create simple system:
ground = mbs.AddObject(ObjectGround())
mbs.AddNode(NodePoint())
body = mbs.AddObject(MassPoint(mass=1, nodeNumber=0))
m0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=ground))
m1 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=body))
mbs.AddObject(CartesianSpringDamper(markerNumbers=[m0,m1], stiffness=[100,100,100]))
mbs.AddLoad(LoadForceVector(markerNumber=m1, loadVector=[10,10,10]))
mbs.Assemble()
simulationSettings = exu.SimulationSettings()
simulationSettings.timeIntegration.endTime = 10
success = mbs.SolveStatic(simulationSettings, storeSolver = True)
exu.Print("success =", success)
exu.Print("iterations = ", mbs.sys['staticSolver'].it)
exu.Print("pos=", mbs.GetObjectOutputBody(body,localPosition=[0,0,0],
variableType=exu.OutputVariableType.Position))
NOTE: this function is directly available in MainSystem (mbs); it should be directly called as mbs.SolveStatic(…). For description of the interface, see the MainSystem Python extensions, Function: SolveStatic
Function: SolveDynamic
function description: solves the dynamic mbs problem using simulationSettings and solver type; see MainSolverImplicitSecondOrder, MainSolverImplicitSecondOrder, for further details of the dynamic solver; this function is also available in exudyn (using exudyn.SolveDynamic(…))
input:
mbs: the MainSystem containing the assembled system; note that mbs may be changed upon several runs of this functionsimulationSettings: specific simulation settings out of exu.SimulationSettings(), as described in SolutionSettings; use options for newton, discontinuous settings, etc., from timeIntegration; therein, implicit second order solvers use settings from generalizedAlpha and explict solvers from explicit; be careful with settings, as the influence accuracy (step size!), convergence and performance (see special Performance and ways to speed up computations)solverType: an exudyn.DynamicSolverType for this run, which takes the place of simulationSettings.timeIntegration.solverType while the solver runs and is set back afterwards; None (default): the solver given by simulationSettings.timeIntegration.solverType (default GeneralizedAlpha)updateInitialValues: if True, the results are written to initial values, such at a consecutive simulation uses the results of this simulation as the initial values of the next simulationstoreSolver: if True, the dynamicSolver object is stored in the mbs.sys dictionary as mbs.sys[‘dynamicSolver’], and simulationSettings are stored as mbs.sys[‘simulationSettings’]showHints: show additional hints, if solver failsshowCausingItems: if linear solver fails, this option helps to identify objects, etc. which are related to a singularity in the linearized system matrixautoAssemble: if True: if mbs.systemIsConsistent=False (system is not assembled), call mbs.Assemble() before solver calls
output: (type: bool) returns True, if successful, False if fails; if storeSolver = True, mbs.sys contains dynamicSolver, which allows to investigate solver problems (see Solver substructures and the items described in MainSolverImplicitSecondOrder)
example:
import exudyn as exu
from exudyn.itemInterface import *
SC = exu.SystemContainer()
mbs = SC.AddSystem()
#create simple system:
ground = mbs.AddObject(ObjectGround())
mbs.AddNode(NodePoint())
body = mbs.AddObject(MassPoint(mass=1, nodeNumber=0))
m0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=ground))
m1 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=body))
mbs.AddObject(CartesianSpringDamper(markerNumbers=[m0,m1], stiffness=[100,100,100]))
mbs.AddLoad(LoadForceVector(markerNumber=m1, loadVector=[10,10,10]))
#
mbs.Assemble()
simulationSettings = exu.SimulationSettings()
simulationSettings.timeIntegration.endTime = 10
success = mbs.SolveDynamic(simulationSettings, storeSolver = True)
exu.Print("success =", success)
exu.Print("iterations = ", mbs.sys['dynamicSolver'].it)
exu.Print("pos=", mbs.GetObjectOutputBody(body,localPosition=[0,0,0],
variableType=exu.OutputVariableType.Position))
NOTE: this function is directly available in MainSystem (mbs); it should be directly called as mbs.SolveDynamic(…). For description of the interface, see the MainSystem Python extensions, Function: SolveDynamic
Function: _SolveDynamic
function description: SolveDynamic with the solver type decided and set in the settings
Function: SolverSuccess
SolverSuccess(solverStructure)
function description: return success (True/False) and error message of solver after SolveSteps(…), SolveSystem(…), SolveDynamic(…) or SolveStatic(…) have been called. May also be set if other higher level functions called e.g. SolveSystem(…)
input:
solverStructure: solver structure, as stored in mbs.sys or as created e.g. by exudyn.MainSolverExplicit()
output: (type: List[bool, str]) returns list [success,errorString], in which success=True or False and in case of no success, information is provided in errorString
example:
#assume MainSystem mbs, exu library and simulationSettings:
try:
mbs.SolveDynamic(simulationSettings)
except exu.ExudynError:
[success, msg] = exu.SolverSuccess(mbs.sys['dynamicSolver'])
exu.Print('success=',success)
exu.Print('error message=',msg)
#alternative:
solver=exu.MainSolverImplicitSecondOrder()
...
[success, msg] = exu.SolverSuccess(solver)
Function: ComputeLinearizedSystem
function description: compute linearized system of equations for ODE2 part of mbs, not considering the effects of algebraic constraints; for computation of eigenvalues and advanced computation with constrained systems, see ComputeODE2Eigenvalues; the current implementation is also able to project into the constrained space, however, this currently does not generally work with non-holonomic systems
input:
mbs: the MainSystem containing the assembled systemsimulationSettings: specific simulation settings used for computation of jacobian (e.g., sparse mode in static solver enables sparse computation)projectIntoConstraintNullspace: if False, algebraic equations (and constraint jacobian) are not considered for the linearized system; if True, the equations are projected into the nullspace of the constraints in the current configuration, using singular value decomposition; in the latter case, the returned list contains [M, K, D, C, N] where C is the constraint jacobian and N is the nullspace matrix (C and N may be an empty list, depending on the following flags)singularValuesTolerance: tolerance used to distinguish between zero and nonzero singular values for algebraic constraints projectionreturnConstraintJacobian: if True, the returned list contains [M, K, D, C, N] where C is the constraint jacobian and N is the nullspace matrix (may be empty)returnConstraintNullspace: if True, the returned list contains [M, K, D, C, N] where C is the constraint jacobian (may be empty) and N is the nullspace matrixautoAssemble: if True: if mbs.systemIsConsistent=False (system is not assembled), call mbs.Assemble() before solver calls
output: (type: [ArrayLike, ArrayLike, ArrayLike]) [M, K, D]; list containing numpy mass matrix M, stiffness matrix K and damping matrix D; for constraints, see options with arguments above, return values may change to [M, K, D, C, N]
notes: consider paper of Agundez, Vallejo, Freire, Mikkola, “The dependent coordinates in the linearization of constrained multibody systems: Handling and elimination”, https://www.sciencedirect.com/science/article/pii/S0020740324000791
example:
import exudyn as exu
from exudyn.utilities import *
import numpy as np
SC = exu.SystemContainer()
mbs = SC.AddSystem()
#
b0 = mbs.CreateMassPoint(referencePosition = [2,0,0],
initialVelocity = [2*0,5,0],
mass = 1, gravity = [0,-9.81,0],
drawSize = 0.5, color=graphics.color.blue)
#
oGround = mbs.AddObject(ObjectGround())
#add vertical spring
oSD = mbs.CreateSpringDamper(itemNumbers=[oGround, b0],
localPosition0=[2,1,0],
localPosition1=[0,0,0],
stiffness=1e4, damping=1e2,
drawSize=0.2)
#
mbs.Assemble()
[M,K,D] = mbs.ComputeLinearizedSystem()
exu.Print('M=\n',M,'\nK=\n',K,'\nD=\n',D)
NOTE: this function is directly available in MainSystem (mbs); it should be directly called as mbs.ComputeLinearizedSystem(…). For description of the interface, see the MainSystem Python extensions, Function: ComputeLinearizedSystem
Function: ComputeODE2Eigenvalues
function description: compute eigenvalues for unconstrained ODE2 part of mbs, which represent the square of the eigenfrequencies (in radiant) of the undamped system; the computation may include constraints in case that ignoreAlgebraicEquations=False (however, this currently does not generally work with non-holonomic systems); for algebraic constraints, however, a dense singular value decomposition of the constraint jacobian is used for the nullspace projection; the computation is done for the initial values of the mbs, independently of previous computations. If you would like to use the current state for the eigenvalue computation, you need to copy the current state to the initial state (using GetSystemState, SetSystemState, see SystemData); note that mass and stiffness matrices are computed in dense mode so far, while eigenvalues are computed according to useSparseSolver.
input:
mbs: the MainSystem containing the assembled systemsimulationSettings: specific simulation settings used for computation of jacobian (e.g., sparse mode in static solver enables sparse computation)useSparseSolver: if False (only for small systems), all eigenvalues are computed in dense mode (slow for large systems!); if True, only the numberOfEigenvalues are computed (numberOfEigenvalues must be set!); Currently, the matrices are exported only in DENSE MODE from mbs, which means that intermediate matrices may become huge for more than 5000 coordinates! NOTE that the sparsesolver accuracy is much less than the dense solvernumberOfEigenvalues: number of eigenvalues and eivenvectors to be computed; if numberOfEigenvalues==0, all eigenvalues will be computed (may be impossible for larger or sparse problems!)constrainedCoordinates: if this list is non-empty (and there are no algebraic equations or ignoreAlgebraicEquations=True), the integer indices represent constrained coordinates of the system, which are fixed during eigenvalue/vector computation; according rows/columns of mass and stiffness matrices are erased; in this case, algebraic equations of the system are ignoredconvert2Frequencies: if True, the square root is computed for eigenvalues, they are converted into frequencies (Hz), and the output is [eigenFrequencies, eigenVectors]useAbsoluteValues: if True, abs(eigenvalues) is used, which avoids problems for small (close to zero) eigenvalues; needed, when converting to frequenciescomputeComplexEigenvalues: if True, the system is converted into a system of first order differential equations, including damping terms; returned eigenvalues are complex and contain the ‘damping’ (=real) part and the eigenfrequency (=complex) part; for this case, set useAbsoluteValues=False (otherwise you will not get the complex values; values are unsorted, however!); also, convert2Frequencies must be False in this case! only implemented for dense solverignoreAlgebraicEquations: if True, algebraic equations (and constraint jacobian) are not considered for eigenvalue computation; otherwise, the solver tries to automatically project the system into the nullspace kernel of the constraint jacobian using a SVD; this gives eigenvalues of the constrained system; eigenvectors are not computedsingularValuesTolerance: tolerance used to distinguish between zero and nonzero singular values for algebraic constraints projectionautoAssemble: if True: if mbs.systemIsConsistent=False (system is not assembled), call mbs.Assemble() before solver calls
output: (type: [ArrayLike, ArrayLike]) [eigenValues, eigenVectors]; eigenValues being a numpy array of eigen values (\(\omega_i^2\), being the squared eigen frequencies in (\(\omega_i\) in rad/s)!), eigenVectors a numpy array containing the eigenvectors in every column
author: Johannes Gerstmayr, Michael Pieber
example:
#take any example from the Examples or TestModels folder, e.g., 'cartesianSpringDamper.py' and run it
#specific example:
import exudyn as exu
from exudyn.utilities import *
import numpy as np
SC = exu.SystemContainer()
mbs = SC.AddSystem()
#
b0 = mbs.CreateMassPoint(referencePosition = [2,0,0],
mass = 1, gravity = [0,-9.81,0],
drawSize = 0.5, color=graphics.color.blue)
#
oGround = mbs.AddObject(ObjectGround())
#add vertical spring
oSD = mbs.CreateSpringDamper(itemNumbers=[oGround, b0],
localPosition0=[2,1,0],
localPosition1=[0,0,0],
stiffness=1e4, damping=1e2,
drawSize=0.2)
#
mbs.Assemble()
#
[eigenvalues, eigenvectors] = mbs.ComputeODE2Eigenvalues()
#==>eigenvalues contain the eigenvalues of the ODE2 part of the system in the current configuration
#
#compute eigenfrequencies in Hz (analytical: 100/2/pi Hz for y-direction):
[eigenvaluesHz, ev] = mbs.ComputeODE2Eigenvalues(convert2Frequencies=True)
#
#compute complex eigenvalues:
[eigenvaluesComplex, ev] = mbs.ComputeODE2Eigenvalues(computeComplexEigenvalues=True,
useAbsoluteValues=False)
NOTE: this function is directly available in MainSystem (mbs); it should be directly called as mbs.ComputeODE2Eigenvalues(…). For description of the interface, see the MainSystem Python extensions, Function: ComputeODE2Eigenvalues
Function: ComputeSystemDegreeOfFreedom
function description: compute system DOF numerically, considering Grübler-Kutzbach formula as well as redundant constraints; uses numpy matrix rank or singular value decomposition of scipy (useSVD=True)
input:
mbs: MainSystem for which DOF shall be computedsimulationSettings: used e.g. for settings regarding numerical differentiation; default settings may be used in most casesthreshold: threshold factor for singular values which estimate the redundant constraintsuseSVD: use singular value decomposition directly, also showing SVD values if verbose=Trueverbose: if True, it will show the singular values and one may decide if the threshold shall be adaptedautoAssemble: if True: if mbs.systemIsConsistent=False (system is not assembled), call mbs.Assemble() before solver calls
output: (type: dict) returns dictionary with key words ‘degreeOfFreedom’, ‘redundantConstraints’, ‘nODE2’, ‘nODE1’, ‘nAE’, ‘nPureAE’, where: degreeOfFreedom = the system degree of freedom computed numerically, redundantConstraints=the number of redundant constraints, nODE2=number of ODE2 coordinates, nODE1=number of ODE1 coordinates, nAE=total number of constraints, nPureAE=number of constraints on algebraic variables (e.g., lambda=0) that are not coupled to ODE2 coordinates
notes: this approach could possibly fail with special constraints! Currently only works with dense matrices, thus it will be slow for larger systems
example:
import exudyn as exu
from exudyn.utilities import *
import numpy as np
SC = exu.SystemContainer()
mbs = SC.AddSystem()
#
b0 = mbs.CreateRigidBody(inertia = InertiaCuboid(density=5000,
sideLengths=[1,0.1,0.1]),
referencePosition = [6,0,0],
initialAngularVelocity = [0,8,0],
gravity = [0,-9.81,0],
graphicsDataList = [exu.graphics.Brick(size=[1,0.1,0.1],
color=graphics.color.orange)])
oGround = mbs.AddObject(ObjectGround())
mbs.CreateGenericJoint(itemNumbers=[oGround, b0], position=[5.5,0,0],
constrainedAxes=[1,1,1, 1,0,0],
rotationMatrixAxes=RotationMatrixX(0.125*pi), #tilt axes
useGlobalFrame=True, axesRadius=0.02, axesLength=0.2)
#
mbs.Assemble()
dof = mbs.ComputeSystemDegreeOfFreedom(verbose=1)['degreeOfFreedom'] #print out details
NOTE: this function is directly available in MainSystem (mbs); it should be directly called as mbs.ComputeSystemDegreeOfFreedom(…). For description of the interface, see the MainSystem Python extensions, Function: ComputeSystemDegreeOfFreedom
Function: CheckSolverInfoStatistics
CheckSolverInfoStatistics(solverName, infoStat, numberOfEvaluations)
function description: helper function for solvers to check e.g. if high number of memory allocations happened during simulation This can happen, if large amount of sensors are attached and output is written in every time step
input: stat=exudyn.special.InfoStat() from previous step, numberOfEvaluations is a counter which is proportional to number of RHS evaluations in method