Written by: Garron Fish – Chief Engineer & Sas Harrison – Systems Engineer
Introduction
This blog entry describes some of the model translation process, focusing on the symbolic processing that is performed.
The symbolic processing of equations is applied to the complete set of equations for the model. This “flat” form of the Modelica model is extracted in the first step of the translation process, an example of a simple electrical model is shown in Figure 1.
Figure 1: Simple electrical circuit with Modelica code extracted
In the translation procedure this flat code is translated into procedural code that serves to facilitate function evaluations to be called as the software works through each time step. This process is summarised below (see Dymola User Manual Volume 2 Section 8.8 and Section 8.9 for a more detailed description of this process). The process however does not limit itself to extracting the equations and putting them in a procedural form.
Assessment of Problem Structure by Dymola
“Structure Jacobian”
From a model’s Modelica flat code (as in Figure 1) a matrix that describes the dependency of each equation on its variables & their derivatives (treated as further variables) is created. The matrix results from each equation being assigned a row and each variable being assigned a column. If an equation makes use of a certain variable then the relevant entry is assigned the value 1 (for example in Figure 2). The only possible values for entries are zero, and unity. This matrix is referred to as the “Structure Jacobian”.
Figure 2: “Structure Jacobian” for a mechanical system. Each row is an equation and each column is a variable.
Image from Dymola User Manual Volume 2.
The “Structure Jacobian” describes the dependency of the equations on the variables; the name “Structure Jacobian” is misleading as it is not a Jacobian, see note at end of post for further details.
The aim is to rearrange the “Structure Jacobian” so that working downwards, each term can be solved for. This rearrangement essentially represents the automated writing of instructions for the computer.
Removing constants and alias reduction
To improve the computational efficiency of the solver, the system of equations is simplified as much as possible. Two simplification methods are applied, with a view to removing constants and “alias” (essentially “duplicate”) variables by substituting at translation.
Firstly, from examining the “Structure Jacobian” it is possible to determine which variables are really constants. This automates the process a person might use to spot that some “variables” for a particular set of boundary conditions or in a particular model configuration, do not vary in time, simplifying the specific solution of the system for that set of parameters & initial conditions. Some software would waste time “calculating” these variables at each time-step. This doesn’t necessarily take long, but it adds up, and blocks up spaces in the memory, which for large systems is very damaging. Dymola avoids it.
Secondly, by examining rows in the “Structure Jacobian” that contain only two variables it is possible to determine which variables may be alias or the negated alias (i.e. b = a or c = -a). The alias variables are removed, for example in Figure 3. This automates the process a person might use to spot quantities that are identically equal, so they don’t end up solving for the same thing twice, needlessly increasing the work. By doing this algorithmically, Dymola lets it happen for potentially hundreds of thousands of equations, so that problems can use the computational speed & efficiency of a computer, while not losing the common-sense of a person.
Figure 3: “Structure Jacobian” example with alias equations removed. The number of unknowns is reduced from about 1200 to 330. Image from Dymola User Manual Volume 2.
A significant reduction in complexity is typically achieved by alias reduction.
DAE to ODE
The Modelica flat code can describe a mixed system of Differential-Algebraic Equations (DAEs – https://en.wikipedia.org/wiki/Differential_algebraic_equation). However DAE solvers are typically slow and not robust with regards to initialisation, so the system is generally transformed into a set of Ordinary Differential Equations (ODEs), to take advantage of a wider variety of numerical solvers.
This transformation is called Index Reduction and a key part of this process is selection of states. The process followed by Dymola supports “dynamic states”, and is described in https://www.modelica.org/events/workshop2000/proceedings/old/Mattsson.pdf. This means that the variables that are considered states in the resulting set of ODEs can be selected dynamically during simulation.
The term “state” is used in state-space modelling, and in principle means “variables chosen to be the coordinates”. For example, let’s imagine you are calculating some dynamics in cartesian space, and then you notice it would be easier or better to use polar coordinates. By changing to doing the equations in a new set of coordinates you are changing the chosen “states” that (fully) describe the configuration of the dynamical system (consider looking at http://web.mit.edu/2.14/www/Handouts/StateSpace.pdf or similar material on state-space modelling).
When Dymola uses static states, Dymola takes a system where the cartesian and polar coordinates are both defined, and chooses which set of variables (polar or cartesian) should be used as the coordinates, and which should just be back-calculated after the fact. This can have a huge effect on the computational efficiency, and is a powerful capability that aids analysis.
The use of dynamic states means that Dymola not only chooses states at the beginning, but continues to re-evaluate whether it should stick to the initial choice, or change, during the analysis! This is similar to how you might solve a motion problem piecemeal, with cartesian coordinates for one part, then switch to polar for another segment which lends itself better to that.
Note that dynamic state selection is typically turned off & requires careful management by a skilled user. Discussion of how these can be used for MultiBody systems can be found here: Understanding State Selection in Multibody Simulation.
Block Lower Triangle Partitioning
In the “ideal” case the Structure Jacobian can be rearranged into a Lower Triangular Matrix. Variables can then be solved for one by one. In general, this is not the case however, and “algebraic loops” (i.e. “coupled equations”) are present in the model. In this case, the Structure Jacobian cannot be restructured as a Lower Triangular Matrix but can be made to approach that state, laid out as in the example in Figure 4.
Figure 4: “structure” Jacobian restructured as Block Lower Triangular Matrix
The blocks along the diagonal of the BLT structure represent algebraic loops. These “algebraic loops” are two-way coupled sets of equations, possibly, but not necessarily, differential equations, requiring simultaneous solution. These could be either of systems of linear equations or systems with some nonlinear equations. The key point to note is that the distinct “blocks” that can be identified, have been isolated in the arrangement of Figure 4.
Separating these blocks, using algebra, already makes the problem more tractable, particularly if one or more of the isolated blocks requires the application of iterative nonlinear solvers. This is because it is much easier if these computationally expensive solvers are only applied to the blocks of equations that need them, and as piecemeal as possible. For example, looking at Figure 4, applying the nonlinear solver to the blocks shown, instead of the whole diagram would be faster. Applying the solver to each block identified in turn, instead of all at once, will be even faster, so the machine doing this algebra ahead of time can really pay off.
Linear system
The linear system of differential equations has the form shown below. Note that the software will reduce higher order systems to this “first order” form, by creating the requisite additional equations in the same way as is documented in the typical elementary undergraduate differential equation courses provided to engineers.
Ax = b
In this equation, the terms of x may well be derivatives of other terms.
An example of a linear system derived from a basic transformer is in Figure 5.
Figure 5: Basic transformer based on Modelica.Electrical.Analog.Examples.CompareTransformers. The alias variables and constants are not shown. The Solve function calculates the solution for der(transformer.i1) and der(transformer.i2) for the given A and b matrices.
In Figure 5 the linear system of equations are restructured into the Ax = b form and x is solved for using the Dymola Solve function.
Nonlinear system of differential equations
A nonlinear system of differential equations has the following form:
fn(x) = 0
where the function fn() has nonlinear characteristics.
This system of equations is solved using a root finding method similar to the Newton-Raphson method (a modified Powell Hybrid method is used). This method iteratively updates x (the iterative variables) based on the gradient of fn(x) until fn(x) converges on 0 (i.e. the residuals converge to zero).
The gradient of fn(x) is the Jacobian matrix (as entirely distinct, conceptually, from the “Structure Jacobian” discussed above):
Figure 6: Jacobian matrix
In the Newton-Raphson method the following update equation is used:
This equation describes how the inverse of the Jacobian matrix is used to update the iterative variables (refer to http://fourier.eng.hmc.edu/e176/lectures/NM/node21.html for more detail), a similar method is used in Dymola.
An example of a nonlinear system of differential equations is in Figure 7. The nonlinear characteristics due to the discontinuity in the saturating inductor.
Figure 7: Model of saturating Inductance driven by sine wave voltage. The Jacobian is calculated symbolically and is used in an iterative process to solve for Inductance.Lact and Inductance.i (the iterative process is not show).
Numerical Jacobians
The iterative process for finding the solution of the nonlinear equations uses the Jacobian matrix of the nonlinear system. Often this Jacobian matrix is calculated symbolically during symbolic processing (as in Figure 7). However sometimes the Jacobian cannot be symbolically calculated by the translation process, in this case the Jacobian matrix is calculated numerically using the forward difference method (see http://www.iue.tuwien.ac.at/phd/khalil/node14.html for more details). This calculation is obviously repeated at each evaluation of the function.
Evaluation of a numerical Jacobian is computationally expensive and the resulting derivative can suffer from noise and inaccuracies.
Depending on the kind of problem & its local behaviour, a numerical derivative can be considerably less satisfactory than a value from a closed-form expression too. If possible, the need for numerical evaluation of derivatives should be avoided (see https://www.claytex.com/blog/how-can-i-make-my-models-run-faster/ and Dymola User Manual Volume 2 Section 8.9 for further details). Take note, that avoiding numerical differentiation while converging nonlinear solutions, can be seen as one key reason to use Dymola software!
Reduction of algebraic loops
Algebraic loops are coupled systems of equations (whether linear or not) which are relatively computationally intensive. There are a number of methods that reduce the number of such coupled systems in the model. Firstly the process of converting the “Structure Jacobian” to the BLT ODE matrix uses efficient graph theory algorithms to select states that generate minimal algebraic loops. This is essentially an automated exploration for appropriate “changes of variable” (see above “Dynamic States”) – you are already familiar with the idea of transformations that decouple equations, and this is the same idea.
Some small linear algebraic loops can be solved symbolically.
Tearing
Tearing is used to reduce the size of algebraic loops. This method can be seen as breaking an algebraic loop by removing the algebraic connection between the output of the system and the input as in Figure 8 below and compensating, by constraining the resulting system such that x2 must equal NEWx2. This is essentially a procedure for iterative solution, particularly likely to be used on coupled systems (called “loops” in the interface, as discussed) which have some nonlinear properties. On coupled sets of equations for which tearing is deployed by the software, iteration variables (called “tearing variables”) are selected intelligently, and a process of converging iterations is undertaken.
Figure 8: Tearing of algebraic loops diagram. Red line is where the f1 and f2 functions are torn
Tearing the model as in Figure 8, allows the remaining equations to be simplified further. See http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.571.9335&rep=rep1&type=pdf for more details.
Other methods for reducing the number of nonlinear systems are (See Section 8.9 of Dymola User Manual Volume 2):
- For scalar nonlinear equations applying inverses functions can sometimes be used to solve the equation analytically.
- Use of general solutions for certain classes of algebraic equations (like polynomials).
- Partitioning of a system of equations into a linear and a nonlinear part can reduce the demands of the nonlinear part.
- Using min and max values to evaluate if-conditions.
Figure 9: Size of algebraic loops is reduced by using tearing
Translation Tab Statistics
Figure 10: A model’s Translation tab in Simulation log of a vehicle model
During the translation process the Simulation log’s translation tab is generated that contains information about symbolic processing that took place.
The Translation log is split into different sections, as in Figure 10, these are:
- A warnings and information section – A number of useful warnings, information and possibly errors.
- A Statistics section – See the next section
- A Settings section – Settings used in the model
- Selected continuous state time states – States selected by the model, or sets from which dynamic states can be selected from
- Warnings regarding initial values of iteration variables – Solutions for nonlinear equations can sometimes be difficult or slow to find. To improve robustness and efficiency of your model you should supply suitable start values for these variables. If no initial value is supplied for a variable, a warning is generated and the initial value is set to zero.
Statistics on the translation log
Figure 11, below, contains the statistics of a vehicle model that is under development.
Figure 11: An example of the translation log statistics expanded out
An explanation of the Statistics in Figure 11 follows:
In the Original Model section the following information is provided:
- Number of components – the number of declared components in the model
- Variables – the number of individual variables in the model
- Constants / Parameters / Unknowns – the number of variables divided by whether they are declared as constant, parameter, or time-varying
- Differentiated variables – the number of variables for which a der(x) statement appears in the model
- Equations / Nontrivial – the number of equations, and the number of this that are “non-trivial” – equations that are “trivial” are usually those generated by connect statements
In the Translated Model section the following information is provided:
- Constants – are constants in the model including those determined during symbolic processing – compare this figure with the number of constants declared in the model
- Free parameters – are parameters that can be altered between simulations without re-translation
- Parameter depending – are parameters dependent on other parameters
- Outputs – outputs of the model
- Continuous time states – states present in the model – this gives an indication of the overall size of model, compare this figure with the number of differentiated variables
- Time-varying variables – variables that vary with time, this includes the number of states
- Alias variables – variables equal to time varying-variables or states, or the negation of these
- Assumed default initial conditions – this is the number of variables for which the “start” value is not fixed
- Number of mixed/discrete systems of equations – the number of linear or nonlinear systems of equations that also include discrete variables
- Sizes of linear systems of equations – the linear systems of equations that are determined by the BLT partitioning
- Sizes after manipulation of the linear systems – the linear system are reduced using the methods summarised in “Reduction of algebraic loops”
- Sizes of nonlinear systems of equations – the nonlinear systems of equations that are determined by the BLT partitioning
- Sizes after manipulation of the nonlinear systems – the nonlinear systems of equations are reduced using the methods summarised in “Reduction of algebraic loops”
- Number of numerical Jacobians – is the number of nonlinear systems that require a numerical Jacobians to be calculated
In Modelica models there are initial equation and initial algorithm sections that can be used to provide characteristics to the initialisation problem, these equations are used together with the system’s equations to characterise the Initialisation system of equations. The Initialisation problem describes the system of initialisation equations and uses the same terminology as used in the Translated model section.
Future blog posts will give further tips on how to understand more about your model using these statistics.
Important Additional Note on the term “Structure Jacobian”
Please take note, that the use of the term “Structure Jacobian” in this blog post, is included to remain consistent with the terminology used by the developers & the manual, however the “Structure Jacobian” concept is not as immediately related to the usual meaning of the term “Jacobian” which will also be used further on, as one might think, and the term is used for reasons beyond the current scope. If it helps, simply read “Structure Array”.
Mathematicians usually use “Jacobian” to mean: The matrix of first-order partial derivatives of a vector function or in other words, a kind of generalised first order derivative that applies to vector functions.
Please get in touch if you have any questions or have got a topic in mind that you would like us to write about. You can submit your questions / topics via: Tech Blog Questions / Topic Suggestion.