Creating multidimensional SD Models

Description and overview how vector-valued and matrix-valued SD Models can bet setup and worked with

Keywords

system dynamics, systemdynamics, sd dsl, bptk, bptk-py, python, business simulation

Multidimensional SD Models

This document is the reference for arrays in the SD DSL: how to give an element a shape, which operators work on it, what each one returns, and what is deliberately not supported. Every section runs its own small example, so you can read it as a catalogue and copy from it.

For arrays in a model that does something, the Model Library has two worked examples: a workforce aging chain over a vector of seniority levels, and a regional product portfolio built on a matrix of products across regions.

We start with some boilerplate to get a BPTK project up and running:

This is already enough to define arrayed components.

How to define arrayed components

There are two options for arrayed components:

• Vectors (one dimensional arrays)
• Matrices (two dimensional arrays)

Moreover, both types of arrays - Vectors and Matrices - can be setup:

• using numerical indices
• using string-valued indices (named arrays)

Lets have a look at some examples:

As can be seen, we need two parameters for setting up a Vector using numerical indices. Moreover, there is one optional parameter.

Parameter	Type	Meaning
size	Integer	Defines the length of the Vector
values	List of Float/Integer	Defines the values of the Vector elements
set_stack_equation	Boolean	(optional) If the element is a stock, the initial value is set (False) or the equation is set (True). Default is False.

And we need one parameter (+ one optional parameter) for setting up a Vector using string indices:

Parameter	Type	Meaning
values	Dictionary	Defines the string-values indices and their values
set_stack_equation	Boolean	(optional) If the element is a stock, the initial value is set (False) or the equation is set (True). Default is False.

values is either one value per index, or a single value that every index gets - setup_vector(2, 3.0) above is the short form of setup_vector(2, [3.0, 3.0]). set_stack_equation has a section of its own below.

For matrices, we can proceed completely similar.

As can be seen, we need two parameters (+ one optional parameter) for setting up a Matrix using numerical indices:

Parameter	Type	Meaning
size	List (tuple) of Integer	Defines the size of the Matrix
values	List of Lists of Float/Integer	Defines the values of the Matrix elements
set_stack_equation	Boolean	(optional) If the element is a stock, the initial value is set (False) or the equation is set (True). Default is False.

And we need one parameter (+ one optional parameter) for setting up a Matrix using string-valued indices:

Parameter	Type	Meaning
values	Dictionary	Defines the string-values indices and their values
set_stack_equation	Boolean	(optional) If the element is a stock, the initial value is set (False) or the equation is set (True). Default is False.

A matrix takes a single value as well: setup_matrix([2, 2], 3.0) gives a 2x2 matrix whose four cells all hold 3.0.

What set_stack_equation does

For every element except a stock, the values handed to setup_vector and its three siblings are simply the sub-elements’ values, and set_stack_equation changes nothing.

A stock is the one case where those values could mean two different things, and this is the flag that decides which:

• False, the default: each value is the initial value of its sub-stock. The sub-stock starts there and then accumulates whatever flows you give it.
• True: each value becomes the sub-stock’s equation, that is its net change per unit of time. The sub-stock starts at 0.0 and adds that value at every step.

So the same call means “start at 10 and 20” or “grow by 10 and 20 per period”:

set_stack_equation=False (initial values, plus a flow of 1.0 and 2.0):
  index 0: [10.0, 11.0, 12.0, 13.0]
  index 1: [20.0, 22.0, 24.0, 26.0]
set_stack_equation=True (the values are the stock's equation):
  index 0: [0.0, 10.0, 20.0, 30.0]
  index 1: [0.0, 20.0, 40.0, 60.0]

Math-Operations for arrayed Components

Arrays are not a special case in the SD DSL: every operator and function that takes an operand works on an arrayed element, element by element. That covers the four arithmetic operations, powers and modulo, the math functions, comparisons and conditionals, max and min, the table and time functions, and the stateful smooth, trend and delay.

On top of that there are the array-specific operations, which are the ones that do something an element-wise operator cannot: they aggregate a whole array into a single value, or multiply arrays in the linear-algebra sense.

Standard Operations

The arithmetic operations \(+\), \(-\), \(*\), \(/\) as well as \(**\) (power) and \(\%\) (modulo) accept any of these operand pairings:

Operand 1	Operand 2
Arrayed Element	Arrayed Element
Arrayed Element	Scalar Element
Scalar Element	Arrayed Element
Arrayed Element	Float/Integer
Float/Integer	Arrayed Element

Scalar Element and Float/Integer are not the same thing. A scalar element is a model element without a shape - a constant, converter, stock or flow - so it has an equation of its own and its value can change over the course of a simulation, and changing it changes every index that reads it. A Float or Integer is a plain Python number written into the equation, fixed for the whole run. Both are allowed on either side of the operator; v * 2.0 and v * some_constant differ only in whether the factor can still move.

⚠️ If both operands are arrayed elements, they must have the same numerical or string-valued indices.

That means it is not possible to have operand 1 = vector with numerical indices and operand 2 = vector with string-valued indices, even if they have the same size.

It is also not possible to have operand 1 = vector and operand 2 = matrix or vice versa - there is no broadcasting. A mismatch raises an exception when the equation is assigned rather than guessing what was meant: mixing a vector and a matrix gives Attempted invalid array addition, two vectors of different length Cannot perform binary operation on arrays with different sizes, and a numerical against a named vector Cannot perform binary operation on arrays with different indices.

Every one of these operations is performed element-wise: index by index, never across indices. Let’s have a look at some examples:

Addition (\(+\))

\[\begin{equation*} \begin{pmatrix} 1.1 \\ 2.2 \end{pmatrix} + \begin{pmatrix} 3.1 \\ 4.2 \end{pmatrix} = \begin{pmatrix} 4.2 \\ 6.4 \end{pmatrix} \end{equation*}\]

[ 4.2 , 6.4 ]

\[\begin{equation*} \begin{pmatrix} 3.1 \\ 4.2 \end{pmatrix} +1.0 = \begin{pmatrix} 4.1 \\ 5.2 \end{pmatrix} \end{equation*}\]

[ 4.1 , 5.2 ]

Subtraction (\(-\))

\[\begin{equation*} \begin{pmatrix} 1.1 & 2.2 \\ 3.3 & 4.4 \end{pmatrix} - \begin{pmatrix} 5.5 & 7.7 \\ 3.3 & 14.4 \end{pmatrix} = \begin{pmatrix} -4.4 & -5.5\\ 0.0 & -10.0 \end{pmatrix} \end{equation*}\]

[ [-4.4 , -5.5]
  [0.0 , -10.0] ]

\[\begin{equation*} \begin{pmatrix} 5.5 & 7.7 \\ 3.3 & 14.4 \end{pmatrix} -1.0 = \begin{pmatrix} 4.5 & 6.7\\ 2.3 & 13.4 \end{pmatrix} \end{equation*}\]

[ [4.5 , 6.7]
  [2.3 , 13.4] ]

Multiplication (\(*\))

\[\begin{equation*} \begin{pmatrix} 4.0\\ 5.0 \end{pmatrix} \odot \begin{pmatrix} 6.0\\ 7.0 \end{pmatrix} = \begin{pmatrix} 24.0\\ 35.0 \end{pmatrix} \end{equation*}\]

[ 24.0 , 35.0 ]

\[\begin{equation*} \begin{pmatrix} 6.0\\ 7.0 \end{pmatrix} \cdot 3.0 = \begin{pmatrix} 18.0\\ 21.0 \end{pmatrix} \end{equation*}\]

[ 18.0 , 21.0 ]

The case “- arrayed element” is a special case since it is interpreted as “(-1) \(\cdot\) element”:

\[\begin{equation*} - \begin{pmatrix} 6.0\\ 7.0 \end{pmatrix} = \begin{pmatrix} -6.0\\ -7.0 \end{pmatrix} \end{equation*}\]

[ -6.0 , -7.0 ]

Division (\(/\))

\[\begin{equation*} \begin{pmatrix} 2.0 & 4.0\\ 8.0 & 16.0 \end{pmatrix} \oslash \begin{pmatrix} 2.0 & 1.0\\ 0.5 & 0.25 \end{pmatrix} = \begin{pmatrix} 1.0 & 4.0\\ 16.0 & 64.0 \end{pmatrix} \end{equation*}\]

[ [1.0 , 4.0]
  [16.0 , 64.0] ]

\[\begin{equation*} \begin{pmatrix} 2.0 & 1.0\\ 0.5 & 0.25 \end{pmatrix} / \text{ } 5.0 = \begin{pmatrix} 0.4 & 0.2\\ 0.1 & 0.05 \end{pmatrix} \end{equation*}\]

[ [0.4 , 0.2]
  [0.1 , 0.05] ]

Power (\(**\)) and Modulo (\(\%\))

This section and the ones that follow all use the same little named vector, so the results are easy to compare.

\[\begin{equation*} \begin{pmatrix} 4.0 \\ 9.0 \end{pmatrix} ^{2} = \begin{pmatrix} 16.0 \\ 81.0 \end{pmatrix} \end{equation*}\]

\[\begin{equation*} \begin{pmatrix} 4.0 \\ 9.0 \end{pmatrix} \bmod 5 = \begin{pmatrix} 4.0 \\ 4.0 \end{pmatrix} \end{equation*}\]

v ** 2: [ 16.0 , 81.0 ]
v % 5:   [ 4.0 , 4.0 ]

Math Functions

Every function in sd_functions that takes a value works on an array and returns an array of the same shape: sqrt, exp, ln, log10, sin, cos, tan, arcsin, arccos, arctan, floor, ceil, round, abs, sinwave, coswave.

\[\begin{equation*} \sqrt{ \begin{pmatrix} 4.0 \\ 9.0 \end{pmatrix} } = \begin{pmatrix} 2.0 \\ 3.0 \end{pmatrix} \end{equation*}\]

\[\begin{equation*} 2 \cdot \ln \begin{pmatrix} 4.0 \\ 9.0 \end{pmatrix} = \begin{pmatrix} 2.7726 \\ 4.3944 \end{pmatrix} \end{equation*}\]

sqrt(v):   [ 2.0 , 3.0 ]
2 * ln(v): [ 2.7726 , 4.3944 ]

Comparisons, If, And, Or, Not

A comparison over an array gives one boolean per index, and If then chooses per index as well. The condition, the then branch and the else branch may each be arrayed or scalar in any combination, as long as the arrayed ones agree on shape.

\[\begin{equation*} \mathrm{If}\left( \begin{pmatrix} 4.0 \\ 9.0 \end{pmatrix} > 6, \begin{pmatrix} 4.0 \\ 9.0 \end{pmatrix} , 0 \right) = \begin{pmatrix} 0.0 \\ 9.0 \end{pmatrix} \end{equation*}\]

If(v > 6, v, 0):     [ 0.0 , 9.0 ]
And(v > 3, v < 6):   [ True , False ]

max and min

sd.max and sd.min compare two operands, element by element. Not to be confused with arr_max and arr_min further down, which take the largest or the smallest value within one array.

\[\begin{equation*} \max\left( \begin{pmatrix} 4.0 \\ 9.0 \end{pmatrix} , 6 \right) = \begin{pmatrix} 6.0 \\ 9.0 \end{pmatrix} \end{equation*}\]

max(v, 6): [ 6.0 , 9.0 ]

smooth, trend and delay

The stateful functions work over arrays as well, and each index carries its own history: smoothing a vector is not one smoothed value copied across the indices, it is one smoothing chain per index. Below, both indices start at the initial value 1.0 and each converges towards its own input - the recurrence \(s_{t+1} = s_t + \frac{dt}{\tau}(x - s_t)\) runs once per index:

\[\begin{equation*} \mathrm{smooth}\left( \begin{pmatrix} 4.0 \\ 9.0 \end{pmatrix} , \tau = 3, s_0 = 1 \right) \Bigr|_{t=4} = \begin{pmatrix} 3.4074 \\ 7.4198 \end{pmatrix} \end{equation*}\]

t= 0: [ 1.0 , 1.0 ]
t= 1: [ 2.0 , 3.6667 ]
t= 4: [ 3.4074 , 7.4198 ]
t=10: [ 3.948 , 8.8613 ]

Array-specific Operations

Array Sum

Calculates the element-wise sum of an array.

\[\begin{equation*} \text{sum} \begin{pmatrix} 1.0 \\ 2.0\\ 3.0 \end{pmatrix} = 1.0 + 2.0 + 3.0 = 6.0 \end{equation*}\]

6.0

Array Product

Calculates the element-wise product of an array.

\[\begin{equation*} \text{prod} \begin{pmatrix} 2.0 & 3.0 & 4.0 \\ 5.0 & 6.0 & 7.0 \end{pmatrix} = 2.0 \cdot 3.0 \cdot 4.0 \cdot 5.0 \cdot 6.0 \cdot 7.0 = 5040.0 \end{equation*}\]

5040.0

Array Rank

Calculates the \(n\)-th highest element of an array. \(n\) is given by a parameter.

For \(n=-1\) the lowest element of the array will be returned.

\[\begin{equation*} \text{rank} \left( \begin{pmatrix} -2.0\\ -0.1\\ 3.1\\ 5.2\\ 11.1 \end{pmatrix} , 1 \right) = 11.1 \end{equation*}\]

11.1

\[\begin{equation*} \text{rank} \left( \begin{pmatrix} -2.0\\ -0.1\\ 3.1\\ 5.2\\ 11.1 \end{pmatrix} , 4 \right) = -0.1 \end{equation*}\]

-0.1

\[\begin{equation*} \text{rank} \left( \begin{pmatrix} -2.0\\ -0.1\\ 3.1\\ 5.2\\ 11.1 \end{pmatrix} , -1 \right) = -2.0 \end{equation*}\]

-2.0

Array Mean

Calculates the element-wise mean of an array.

\[\begin{equation*} \text{mean} \begin{pmatrix} 2.0 & 4.0\\ 6.0 & 8.0\\ 10.0 & 12.0 \end{pmatrix} = \frac{2.0 + 4.0 + 6.0 + 8.0 + 10.0 +12.0}{6} = 7.0 \end{equation*}\]

7.0

Array Median

Calculates the element-wise median of an array.

\[\begin{equation*} \text{median} \begin{pmatrix} -2.0\\ -0.1\\ 3.1\\ 5.2\\ 11.1 \end{pmatrix} = 3.1 \end{equation*}\]

3.1

\[\begin{equation*} \text{median} \begin{pmatrix} -2.0\\ -0.1\\ 3.1\\ 5.2\\ \end{pmatrix} = \frac{-0.1+3.1}{2}=1.5 \end{equation*}\]

1.5

Array Standard Deviation

Calculates the element-wise standard deviation of an array.

\[\begin{equation*} \sigma \begin{pmatrix} 1.0 & 3.0\\ 3.0 & 1.0 \end{pmatrix} = \sqrt{ \frac{1}{4} \cdot \left( (1-2)^2 + (3-2)^2 + (3-2)^2 + (1-2)^2 \right) } = 1 \end{equation*}\]

1.0

Array Maximum and Minimum

arr_max returns the largest value of an array, arr_min the smallest - over every element of a vector or every cell of a matrix.

Both take no arguments, and both return a single value however arrayed the input is. On an element with no sub-elements they return 0.0, as the other aggregations do.

\[\begin{equation*} \mathrm{arr\_max} \begin{pmatrix} 3.0 & 8.0 \\ 1.0 & 5.0 \end{pmatrix} = 8.0 \qquad \mathrm{arr\_min} \begin{pmatrix} 3.0 & 8.0 \\ 1.0 & 5.0 \end{pmatrix} = 1.0 \end{equation*}\]

arr_max: 8.0
arr_min: 1.0

Array Size

Calculates the size of an array.

For a vector, the length will be returned. For a matrix, the size of the highest level will be returned (for example 2 for a \(2 \times 3\) matrix).

\[\begin{equation*} \text{len} \begin{pmatrix} 1.0\\ 1.0\\ 1.0\\ 1.0\\ 1.0\\ 1.0\\ \end{pmatrix} = 6 \end{equation*}\]

6

\[\begin{equation*} \text{len} \begin{pmatrix} 1.0 & 1.0 & 1.0\\ 1.0 & 1.0 & 1.0 \end{pmatrix} = 2 \end{equation*}\]

2

\[\begin{equation*} \text{len} \begin{pmatrix} 1.0 & 1.0 & 1.0\\ 1.0 & 1.0 & 1.0\\ 1.0 & 1.0 & 1.0\\ 1.0 & 1.0 & 1.0 \end{pmatrix} = 4 \end{equation*}\]

4

Array Dot

The Dot function provides the classical vector/matrix-multiplication logic. That means, the following can be calculated:

Factor 1	Factor 2	Result
Vector of size \(m\)	Constant	Vector of size \(m\)
Constant	Vector of size \(m\)	Vector of size \(m\)
Matrix of size \(m \times n\)	Constant	Matrix of size \(m \times n\)
Constant	Matrix of size \(m \times n\)	Matrix of size \(m \times n\)
Vector of size \(m\)	Vector of size \(m\)	Value (Scalar Product)
Vector of size \(m\)	Matrix of size \(m \times n\)	Vector of size \(n\)
Matrix of size \(m \times n\)	Vector of size \(n\)	Vector of size \(m\)
Matrix of size \(m \times n\)	Matrix of size \(n \times p\)	Matrix of size \(m \times p\)

❗ Using the Dot function for an array and a constant yields the same result as using the \(*\)-Operator for the array and the value of the constant.

If the dimensions of the arrays to which the dot function is applied do not allow for a valid array multiplication, an exception is thrown.

⚠️ The Dot function is currently supported for not-named arrays only!

Lets have a look at some examples:

\[\begin{equation*} 2.0 \cdot \begin{pmatrix} 1.0 \\ 2.0 \\ 3.0 \end{pmatrix} = \begin{pmatrix} 2.0 \\ 4.0 \\ 6.0 \end{pmatrix}\end{equation*}\]

[2.0 , 4.0 , 6.0]

\[\begin{equation*} \begin{pmatrix} 4.0 \\ 5.0 \\ 6.0 \end{pmatrix} \cdot 2.0 = \begin{pmatrix} 8.0 \\ 10.0 \\ 12.0 \end{pmatrix}\end{equation*}\]

[8.0 , 10.0 , 12.0]

\[\begin{equation*} \begin{pmatrix} 1.0 & 2.0 \\ 3.0 & 4.0 \\ 5.0 & 6.0 \end{pmatrix} \cdot 2.0 = \begin{pmatrix} 2.0 & 4.0 \\ 6.0 & 8.0 \\ 10.0 & 12.0 \end{pmatrix}\end{equation*}\]

[ [2.0 , 4.0]
  [6.0 , 8.0]
  [10.0 , 12.0] ]

\[\begin{equation*} 2.0 \cdot \begin{pmatrix} -1.0 & -2.0 & -3.0 \\ -4.0 & -5.0 & -6.0 \\ \end{pmatrix} = \begin{pmatrix} -2.0 & -4.0 & -6.0 \\ -8.0 & -10.0 & -12.0 \\ \end{pmatrix}\end{equation*}\]

[ [-2.0 , -4.0 , -6.0]
  [-8.0 , -10.0 , -12.0] ]

\[\begin{equation*} \left\langle \begin{pmatrix} 1.0 \\ 2.0 \\ 3.0 \end{pmatrix}, \begin{pmatrix} 4.0 \\ 5.0 \\ 6.0 \end{pmatrix} \right\rangle = 1.0 \cdot 4.0 + 2.0 \cdot 5.0 + 3.0 \cdot 6.0 = 32.0 \end{equation*}\]

32.0

\[\begin{equation*} \begin{pmatrix} 1.0 & 2.0 & 3.0 \\ \end{pmatrix} \cdot \begin{pmatrix} 1.0 & 2.0 \\ 3.0 & 4.0 \\ 5.0 & 6.0 \\ \end{pmatrix} = \begin{pmatrix} 1.0 \cdot 1.0 + 2.0 \cdot 3.0 + 3.0 \cdot 5.0 \\ 1.0 \cdot 2.0 + 2.0 \cdot 4.0 + 3.0 \cdot 6.0 \\ \end{pmatrix} = \begin{pmatrix} 22.0 & 28.0 \end{pmatrix} \end{equation*}\]

[22.0 , 28.0]

\[\begin{equation*} \begin{pmatrix} -1.0 & -2.0 & -3.0 \\ -4.0 & -5.0 & -6.0 \\ \end{pmatrix} \cdot \begin{pmatrix} 4.0 \\ 5.0 \\ 6.0 \\ \end{pmatrix} = \begin{pmatrix} -1.0 \cdot 4.0 + (-2.0) \cdot 5.0 + (-3.0) \cdot 6.0 \\ -4.0 \cdot 4.0 + (-5.0) \cdot 5.0 + (-6.0) \cdot 6.0 \\ \end{pmatrix} = \begin{pmatrix} -32.0 \\ -77.0 \end{pmatrix} \end{equation*}\]

[-32.0 , -77.0]

\[\begin{equation*} \begin{aligned} \begin{pmatrix} 1.0 & 2.0 \\ 3.0 & 4.0 \\ \end{pmatrix} \cdot \begin{pmatrix} -1.0 & -2.0 \\ -4.0 & -5.0 \\ \end{pmatrix} &= \begin{pmatrix} 1.0 \cdot (-1.0) + 2.0 \cdot (-4.0) & 1.0 \cdot (-2.0) + 2.0 \cdot (-5.0) \\ 3.0 \cdot (-1.0) + 4.0 \cdot (-4.0) & 3.0 \cdot (-2.0) + 4.0 \cdot (-5.0) \\ \end{pmatrix}\\ &= \begin{pmatrix} -9 & -12 \\ -19 & -26\\ \end{pmatrix} \end{aligned} \end{equation*}\]

[ [-9.0 , -12.0]
  [-19.0 , -26.0] ]

Plotting arrayed Components

Similar to one-dimensional SD DSL elements, we can also plot these elements. Lets have a look:

As can be seen, both elements of the Vector are plotted.

If you want to plot the values of the Matrix, you need to specify the first index.

A simple Example

Lets have a look on a concrete example, how a multidimensional SD Model can look like.

Consider an investment depot with two accounts:

• bank account
• depot account

Both accounts will have different deposit rates and different interest rates each year. We want to investigate the value development of the bank account, the depot account and the whole investment depot.

Lets set up the model:

As always we define a scenario manager and scenarios:

And plot the results:

As always we can compare different scenarios with each other by plotting them simultaneously:

What Is Not Supported

Worth knowing before you build on arrays. None of these fails silently - each one raises with a message that says what happened.

• Two dimensions is the limit. setup_matrix takes exactly two sizes, and there is no three-dimensional array.
• dot works on unnamed arrays only. A named operand raises rather than guessing how to line up the labels.
• There is no aggregation over a single dimension. arr_sum and arr_prod take a dimensions argument, but the only values it accepts are "*" (the default) and an integer equal to the array’s depth - both of which aggregate every cell. Aggregating along the rows of a matrix would have to return a vector, which is not supported; where you need it, address the rows yourself. A row of a named matrix is a named vector in its own right, so matrix['north'].arr_sum() is the total of that row.
• arr_size counts the first dimension, not the number of cells: 2 for a \(2 \times 3\) matrix.
• There is no broadcasting. Both operands of an element-wise operation must have the same shape and the same indices; a vector and a matrix is an error.
• No transpose, and no dimension-position operator. The XMILE standard has both - a transpose, and @ to name a position within a dimension - and the SD DSL has never implemented either, because nothing in the model library or the test corpus has needed one. Where you would reach for a transpose, address the cells the other way round when you set the matrix up; a row of a named matrix is a named vector, so matrix['north'] is that row and there is no need to turn the matrix around to get at it.
• A flow never goes negative, arrayed or not: every flow equation is wrapped in max(0, ...). A quantity that has to move both ways belongs in a biflow, which takes the same setup_vector and setup_named_vector as any other element.