SDXorg / pysd

System Dynamics Modeling in Python
http://pysd.readthedocs.org/
MIT License
350 stars 89 forks source link

Definition of the DelayN function #391

Closed universmile closed 1 year ago

universmile commented 1 year ago

The definition of the DELAY 3I function of Vensim is provided in their documentation, including a complete pseudo-code. In contrast, we could not find the precise definition of the DELAY N function. We are interested particularly in the setting where the delay time is variable.

We recently discovered your nice library about system dynamics. The Delay N function in PySD is implemented as a class and it seems that the desired output is produced by the ddt method, but we did not manage to understand precisely how the final output is produced (e.g., with an order of 10, a variable delay time and a constant initial value).

Could someone share the precise mathematical definition of the Delay N function?

enekomartinmartinez commented 1 year ago

Hi @universmile

First, you must know that the DelayN we implemented is based on the test results. Therefore, we cannot assure you that this will work exactly the same as Vensim always, as we may end up with an algorithm that may not be exactly the same but behave the same in our current tests.

The main difference in Vensim between DELAY 3/3I/1/1I and DELAY N depends only on the delay time of the previous time steps. Therefore, when computing the delay, we need to keep a vector of times to ensure the correct delay time is used for each iteration.

Then the implementation of DELAY N is the same that the one used for the others, only that the correspondent delay time is used in each level (if the delay time is constant, there will be no difference).

So the pseudocode for a DELAY N of order 3 could be (if I am not wrong):

DELAY N (order 3)=LV3/DL3
LV3=INTEG(RT2-DELAY3I,
        input*DL1)
DL3=delay time3/3
RT2=LV2/DL1
LV2=INTEG(RT1-RT2,LV3)
DL2=delay time2/3
RT1=LV1/DL1
LV1=INTEG(input-RT1,LV3)
DL1=delay time1/3

why for the regular DELAY 3 (or 3I,1,1I) is:

DELAY 3=LV3/DL3
LV3=INTEG(RT2-DELAY3I,
        input*DL3)
DL3=delay time3/3
RT2=LV2/DL2
LV2=INTEG(RT1-RT2,LV3)
DL2=delay time2/3
RT1=LV1/DL1
LV1=INTEG(input-RT1,LV3)
DL1=delay time1/3

In our implementation, you can compare both; note that the ddt return is the differential to integrate (right argument of INTEG). I hope this helps.

universmile commented 1 year ago

Thank you @enekomartinmartinez for your detailed answer!