diff --git a/docs/whats_new.rst b/docs/whats_new.rst
index 9932c02b..a58eecaf 100644
--- a/docs/whats_new.rst
+++ b/docs/whats_new.rst
@@ -1,7 +1,7 @@
 What's New
 ==========
 
-v3.9.1 (2023/02/xx)
+v3.9.1 (2023/03/11)
 -------------------
 
 New Features
@@ -19,6 +19,7 @@ Bug fixes
 - Set :py:mod:`numpy` <1.24 to avoid errors with least squares equation in :py:func:`pysd.py_backend.allocation.allocate_available`. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
 - Keep the attributes of a component when using :py:meth:`pysd.py_backend.model.Macro.set_components` to avoid losing coords or arguments information. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
 - Set :py:mod:`openpyxl` <3.1 to avoid errors due to non-backwards compatible changes. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+- Include time dependency in random functions to avoid them using constant cache. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
 
 Documentation
 ~~~~~~~~~~~~~
diff --git a/pysd/builders/python/python_expressions_builder.py b/pysd/builders/python/python_expressions_builder.py
index 68bb7db4..64fd8943 100644
--- a/pysd/builders/python/python_expressions_builder.py
+++ b/pysd/builders/python/python_expressions_builder.py
@@ -642,6 +642,7 @@ def build_function_call(self, arguments: dict) -> BuildAST:
                     self.def_subs)[1]
                 expression = f"xr.DataArray({expression}, {subs}, "\
                     f"{list(self.def_subs)})"
+            calls["time"] = 1
 
         elif self.function == "active_initial":
             # Ee need to ensure that active initial outputs are always the
diff --git a/tests/pytest_pysd/pytest_random.py b/tests/pytest_pysd/pytest_random.py
index 7e6e1b1d..37596405 100644
--- a/tests/pytest_pysd/pytest_random.py
+++ b/tests/pytest_pysd/pytest_random.py
@@ -2,6 +2,9 @@
 import pytest
 import shutil
 
+import numpy as np
+import xarray as xr
+
 import pysd
 
 
@@ -28,4 +31,28 @@ def test_translate(self, model_path):
         """
         # expected file
         model = pysd.read_vensim(model_path)
-        model.run()
+        random_vars = [
+            "A B uniform matrix",
+            "A B uniform matrix 1",
+            "A B uniform matrix 1 0",
+            "A B uniform scalar",
+            "A B uniform vec",
+            "A B uniform vec 1",
+            "normal A B uniform matrix",
+            "normal A B uniform matrix 1",
+            "normal A B uniform matrix 1 0",
+            "normal scalar",
+            "normal vec",
+            "normal vec 1",
+            "uniform matrix",
+            "uniform scalar",
+            "uniform vec"
+        ]
+        out = model.run(return_columns=random_vars, flatten_output=False)
+        for var in out.columns:
+            if isinstance(out[var].values[0], xr.DataArray):
+                values = np.array([a.values for a in out[var].values])
+            else:
+                values = out[var].values
+            # assert all values are different in each dimension and time step
+            assert len(np.unique(values)) == np.prod(values.shape)