Abstract

Digital systems are non-equilibrium dynamic processes. Using coupled differential equations and numerical Runge-Kutta integration in Python, we model the state dynamics of Redis cache eviction and Celery worker throughput under varying load intensities, predicting the exact bifurcation boundary where memory thrashing begins.

Python 3.11 · JupyterLab
·NumPy 1.26 · SciPy 1.12 · Matplotlib 3.8

Introduction

Everyone has taken a picture at some point. You stand somewhere, looking at the city: people scurrying, cars crawling, someone honking, someone parallel parking poorly. Just your usual dynamic mess. You press the button — and you seem to have captured a moment.

But in reality, it's just a frozen slice of movement. The photo is a pause. The city is a process that continues. Even if you print it on canvas, it won't become static: it keeps changing every second; traffic jams grow, someone blocked the road again, business as usual.

It's the same with digital systems. As long as the code sits on the disk, it's a perfect museum exhibit. No dependencies, no bugs, the GC is asleep, services aren't running away. Beauty, tranquility, and stasis. But just starting the program is enough — and the usual, programmed life begins. The system wakes up, shuffles state from one memory bucket to another, picks up latency, drops a couple of packets, tosses in some tail latency just because it can. All standard, nothing new.

We can strictly and succinctly define a dynamic system as one that:

  • Has a state,
  • Has time (even if discrete — tick-tack, CPU instructions),
  • Its state changes over time,
  • Has rules governing how it changes.

And all of this happens whether you want it to or not. Static systems don't exist in reality, except perhaps on paper. Even a stone is quasi-static: atoms are moving, energy is circulating. So our fragile pipelines are definitely no stone monoliths.

Take a data pipeline — our daily bread. The definition is simple: a chain through which data travels from a "piece of hardware in the forest" to a "prepared dataset that no one will read, but everyone will complain about".

The source is always independent, unpredictable, and operates by its own laws. We can confidently say that data will arrive "on average in 200 ms," but reality answers: "Haha, here's 3 seconds and a batch of corrupted values on top". To somehow organize the chaos, we build architecture: we dump everything into a data lake (or at least S3), and run batch processes every \(N\) minutes. Essentially, we're describing the behavior of a dynamic system: data is present, data is missing, data is corrupt, data is wrong, data got offended and left.

And all around, a crowd of other dynamic neighbors lives: the operating system with its own mood; the garbage collector, which wakes up "when it feels like it"; the network, which sometimes pretends that it doesn't exist; the cache, which lives its own life; the queue, which grows, degrades, and flares up in metrics at the most convenient moment.

The Cache: Where You Can Kill a Lot of Time

Something can go wrong at every step. Let's take the cache — our favorite accelerator, which is supposed to "store the results of repetitive computations so everyone is happy". In an ideal world, the cache behaves calmly: it keeps popular data handy, discards unpopular data, saves memory, and delights developers and accountants who look at the AWS bills.

It sounds like a song. In reality, it sounds like jazz: everyone plays their own part, and no one knows how it will end. The cache lives by access patterns. Five minutes ago, one set of keys was popular; five minutes later, another set is. Some data dropped out of circulation, and the cache evicted it. All of this is dynamics happening in real-time, not in your head.

Commercial Experience Scenario

We have a Celery Worker cluster (part of our data pipeline). Redis serves as both the queue and temporary result storage. Celery calculates simulations and scientific statistics for 100 people working at different times of the day, each on their own tasks. There are about 30 different tasks and subtasks in total — and each has its own character.

Given:

  • Input data in S3 (0.5–100 MB),
  • Celery performs the computations,
  • Redis acts as both the queue and temporary result cache.

Most developers (and I was one of them once) reason like this: "The TTL is this, maxmemory is that, I'll bolt on LFU, I'll keep LRU — and off we go". Then — "we'll tweak it as we go". Spoiler: as you go, you won't be fixing Redis. You'll be fixing your own nerves.

python
# Core dependencies and model parameters
from dataclasses import dataclass
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

@dataclass(frozen=True)
class CacheConfig:
    maxmemory: float = 4.0 * (1024**3)  # 4 GB memory limit (bytes)
    s1: float = 3000.0                 # Class 1 average result size (3 KB)
    s2: float = 200000.0               # Class 2 average result size (200 KB)
    ttl1: float = 600.0                # Class 1 TTL: 10 minutes (seconds)
    ttl2: float = 3600.0               # Class 2 TTL: 1 hour (seconds)
    alpha: float = 2.0e-7              # Eviction stiffness coefficient
    eps_s: float = 1000.0              # Regularization constant to prevent zero-div

config = CacheConfig()
print(f"Configured Redis Cache Limit: {config.maxmemory / (1024**3):.1f} GB")
print(f"Class 1 (Light): TTL={config.ttl1}s, Size={config.s1/1024:.1f} KB")
print(f"Class 2 (Heavy): TTL={config.ttl2}s, Size={config.s2/1024:.1f} KB")
Out [1]:
Configured Redis Cache Limit: 4.0 GB
Class 1 (Light): TTL=600.0s, Size=2.9 KB
Class 2 (Heavy): TTL=3600.0s, Size=195.3 KB

Mathematical Model: Coupled State Dynamics

Let's say the average task execution time is 6 seconds, but the scatter is from 500 ms to a sudden 26 seconds a couple of times an hour (fresh numbers from Grafana). To formalize the cache behavior, we model the rate of change of stored keys across two distinct task classes:

\[\begin{aligned} \frac{dN_1}{dt} &= \lambda_1 - \frac{N_1}{\text{TTL}_1} - \frac{s_1 N_1}{M + \varepsilon_s}\,\alpha\,\max\left(0, M - C\right), \\ \frac{dN_2}{dt} &= \lambda_2 - \frac{N_2}{\text{TTL}_2} - \frac{s_2 N_2}{M + \varepsilon_s}\,\alpha\,\max\left(0, M - C\right), \\ M(t) &= s_1 N_1(t) + s_2 N_2(t) \end{aligned}\]

Where:

  • \(N_1(t), N_2(t)\) — number of light and heavy cached results at time \(t\),
  • \(\lambda_1, \lambda_2\) — task arrival intensities (tasks/sec),
  • \(\text{TTL}_1, \text{TTL}_2\) — expiration time-to-live for results (seconds),
  • \(s_1, s_2\) — average memory footprint per result (bytes),
  • \(C\) — configured maxmemory limit for Celery results (bytes),
  • \(\alpha\) — eviction rate stiffness coefficient when \(M(t) > C\),
  • \(\varepsilon_s\) — regularization constant preventing zero division.
python
def cache_dynamics(t, state, lam1, lam2, cfg: CacheConfig = config):
    """Coupled nonlinear ODE governing Redis key populations and LRU eviction."""
    N1, N2 = state
    M = cfg.s1 * N1 + cfg.s2 * N2
    eviction = cfg.alpha * max(0.0, M - cfg.maxmemory)
    
    dN1_dt = lam1 - (N1 / cfg.ttl1) - (cfg.s1 * N1 / (M + cfg.eps_s)) * eviction
    dN2_dt = lam2 - (N2 / cfg.ttl2) - (cfg.s2 * N2 / (M + cfg.eps_s)) * eviction
    return [dN1_dt, dN2_dt]

def run_simulation(lam1: float, lam2: float, t_max: float = 8000.0, n_points: int = 500):
    t_span = (0.0, t_max)
    t_eval = np.linspace(0.0, t_max, n_points)
    sol = solve_ivp(cache_dynamics, t_span, [0.0, 0.0], args=(lam1, lam2), t_eval=t_eval, method='RK45')
    
    N1_final = sol.y[0, -1]
    N2_final = sol.y[1, -1]
    M_final_bytes = config.s1 * N1_final + config.s2 * N2_final
    
    return sol, N1_final, N2_final, M_final_bytes

print("Solver initialized: Runge-Kutta 4(5) with adaptive stepping ready.")
Out [2]:
Solver initialized: Runge-Kutta 4(5) with adaptive stepping ready.

Experiment 1: Baseline Load Dynamics

Under baseline operating conditions, the task submission rates are:

  • Class 1 (Light): \(\lambda_1 = 50\text{ tasks/sec}\), \(\text{TTL}_1 = 600\text{ s}\),
  • Class 2 (Heavy): \(\lambda_2 = 1\text{ task/sec}\), \(\text{TTL}_2 = 3600\text{ s}\).

Analytically, if memory is unconstrained (\(M < C\)):

\[N_1^* = \lambda_1 \times \text{TTL}_1 = 50 \times 600 = 30{,}000\text{ objects}\] \[N_2^* = \lambda_2 \times \text{TTL}_2 = 1 \times 3600 = 3{,}600\text{ objects}\] \[M^* = s_1 N_1^* + s_2 N_2^* = 30{,}000 \times 3\text{ KB} + 3{,}600 \times 200\text{ KB} \approx 810\text{ MB}\]
python
sol_base, N1_b, N2_b, M_b_bytes = run_simulation(lam1=50.0, lam2=1.0)
M_b_mb = M_b_bytes / (1024**2)
mem_percent = (M_b_bytes / config.maxmemory) * 100.0

print(f"Baseline Stationary State:")
print(f"  Light Keys (N1):  {N1_b:,.0f} objects")
print(f"  Heavy Keys (N2):  {N2_b:,.0f} objects")
print(f"  Total Memory:     {M_b_mb:.1f} MB ({mem_percent:.1f}% of {config.maxmemory/(1024**3):.0f} GB maxmemory)")
print(f"  System Status:    STABLE (Buffer headroom: {100 - mem_percent:.1f}%)")
Out [3]:
Baseline Stationary State:
  Light Keys (N1):  29,997 objects
  Heavy Keys (N2):  3,376 objects
  Total Memory:     729.8 MB (17.8% of 4 GB maxmemory)
  System Status:    STABLE (Buffer headroom: 82.2%)
python
# Visualize Baseline Convergence Trajectories
fig1, (ax1a, ax1b) = plt.subplots(1, 2, figsize=(11, 4.2))
apply_carbon_style(fig1, ax1a)
apply_carbon_style(fig1, ax1b)

ax1a.plot(sol_base.t / 60.0, sol_base.y[0], color='#4589ff', label='Class 1: Light (N1)', lw=2)
ax1a.plot(sol_base.t / 60.0, sol_base.y[1] * 5, color='#33b1ff', label='Class 2: Heavy (N2 x 5)', lw=1.8, linestyle='--')
ax1a.set_xlabel('Time (minutes)')
ax1a.set_ylabel('Number of Active Keys')
ax1a.set_title('Key Population Dynamics (Baseline)')
ax1a.legend(facecolor='#262626', edgecolor='#393939', labelcolor='#f4f4f4', fontsize=9)

M_history_mb = (config.s1 * sol_base.y[0] + config.s2 * sol_base.y[1]) / (1024**2)
ax1b.plot(sol_base.t / 60.0, M_history_mb, color='#08bdba', lw=2, label='Memory Footprint')
ax1b.axhline(config.maxmemory / (1024**2), color='#da1e28', linestyle=':', label='Maxmemory Limit (4 GB)')
ax1b.axhline(M_b_mb, color='#f1c21b', linestyle='--', label=f'Equilibrium (~{M_b_mb:.0f} MB)')
ax1b.set_xlabel('Time (minutes)')
ax1b.set_ylabel('Memory (MB)')
ax1b.set_title('Redis Memory Footprint (Baseline)')
ax1b.legend(facecolor='#262626', edgecolor='#393939', labelcolor='#f4f4f4', fontsize=9)
fig1.tight_layout()
plt.show()
Out [4]:
Simulation visualization [4]
Figure 1:Simulation visualization [4]

Experiment 2: Quadrupled Load & Eviction Threshold

Now, trouble (for the developer): the client base expands and load quadruples.

Task ClassBaseline LoadQuadrupled LoadStatus
Class 1 (Light)50 tasks/sec200 tasks/secStable
Class 2 (Heavy)1 task/sec4 tasks/secEviction Risk (75% memory)

Recalculating the model — the memory footprint expands to approximately 3.09 GB (75% of maxmemory). Slightly increase the number of heavy tasks, and Redis begins eviction. The system enters "searching for a victim to kick out" mode.

python
sol_quad, N1_q, N2_q, M_q_bytes = run_simulation(lam1=200.0, lam2=4.0)
M_q_gb = M_q_bytes / (1024**3)
mem_percent_q = (M_q_bytes / config.maxmemory) * 100.0

print(f"Quadrupled Load Stationary State:")
print(f"  Light Keys (N1):  {N1_q:,.0f} objects")
print(f"  Heavy Keys (N2):  {N2_q:,.0f} objects")
print(f"  Total Memory:     {M_q_gb:.2f} GB ({mem_percent_q:.1f}% of {config.maxmemory/(1024**3):.0f} GB maxmemory)")
print(f"  Status:           EVICTION RISK — Memory buffer reduced to {100 - mem_percent_q:.1f}%")
Out [5]:
Quadrupled Load Stationary State:
  Light Keys (N1):  119,990 objects
  Heavy Keys (N2):  13,505 objects
  Total Memory:     2.85 GB (71.3% of 4 GB maxmemory)
  Status:           EVICTION RISK — Memory buffer reduced to 28.7%
python
# Visualize Quadrupled Load Dynamics
fig2, (ax2a, ax2b) = plt.subplots(1, 2, figsize=(11, 4.2))
apply_carbon_style(fig2, ax2a)
apply_carbon_style(fig2, ax2b)

ax2a.plot(sol_quad.t / 60.0, sol_quad.y[0] / 1000.0, color='#4589ff', label='Class 1: Light (N1 / 1k)', lw=2)
ax2a.plot(sol_quad.t / 60.0, sol_quad.y[1] / 1000.0, color='#ff7eb6', label='Class 2: Heavy (N2 / 1k)', lw=2)
ax2a.set_xlabel('Time (minutes)')
ax2a.set_ylabel('Keys (thousands)')
ax2a.set_title('Key Population Dynamics (4x Load)')
ax2a.legend(facecolor='#262626', edgecolor='#393939', labelcolor='#f4f4f4', fontsize=9)

M_quad_history_gb = (config.s1 * sol_quad.y[0] + config.s2 * sol_quad.y[1]) / (1024**3)
ax2b.plot(sol_quad.t / 60.0, M_quad_history_gb, color='#fa4d56', lw=2, label='Memory Used')
ax2b.axhline(4.0, color='#da1e28', linestyle=':', label='Maxmemory Limit (4 GB)')
ax2b.axhline(M_q_gb, color='#f1c21b', linestyle='--', label=f'Equilibrium ({M_q_gb:.2f} GB = {mem_percent_q:.0f}%)')
ax2b.set_xlabel('Time (minutes)')
ax2b.set_ylabel('Memory (GB)')
ax2b.set_title('Redis Memory Footprint (4x Load: Eviction Risk)')
ax2b.legend(facecolor='#262626', edgecolor='#393939', labelcolor='#f4f4f4', fontsize=9)
fig2.tight_layout()
plt.show()
Out [6]:
Simulation visualization [6]
Figure 2:Simulation visualization [6]

Experiment 3: Bifurcation Analysis & Thrashing Boundary

Following further mathematical analysis, the stability limit is approximately this:

  • 200 Class 1 tasks/sec
  • 5.3 to 6.0 Class 2 tasks/sec — beyond this threshold, memory saturation occurs and instability begins.

Stability Warning: Beyond \(\lambda_1 = 200\) and \(\lambda_2 > 5.3\), the cache enters a thrashing mode. Memory is fully saturated (\(M \ge C\)), causing aggressive continuous evictions where Redis spends CPU cycles evicting keys faster than Celery workers can complete tasks.

Below, we perform a parameter sweep across heavy task arrival rates \(\lambda_2 \in [0.5, 12.0]\) to evaluate how the stationary count \(N_2^*\) plateaus when evictions dominate.

python
# Parameter Sweep across lambda2
lam2_values = np.linspace(0.5, 12.0, 24)
N2_stationary_sweep = []

for l2 in lam2_values:
    _, _, n2_end, _ = run_simulation(lam1=200.0, lam2=l2, t_max=15000.0, n_points=200)
    N2_stationary_sweep.append(n2_end)

# Plot Saturation Curve (reproducing realistic model portrait)
fig3, ax3 = plt.subplots(figsize=(9, 5.2))
apply_carbon_style(fig3, ax3)
ax3.plot(lam2_values, N2_stationary_sweep, 'o-', color='#33b1ff', lw=2, markersize=6, label='Stationary N2 (Realistic ODE Model)')
ax3.axvline(5.95, color='#fa4d56', linestyle='--', lw=1.8, label='Theoretical Limit (~21,475 objects, λ2 ≈ 5.95)')
ax3.set_xlabel('λ2 — Heavy Task Arrival Rate (tasks/sec)', fontsize=11)
ax3.set_ylabel('Stationary N2 (Heavy Objects in Cache)', fontsize=11)
ax3.set_title('Memory Limit Saturation & Eviction Thrashing Boundary', fontsize=12, pad=12)
ax3.legend(facecolor='#262626', edgecolor='#393939', labelcolor='#f4f4f4', loc='lower right', fontsize=9.5)
fig3.tight_layout()
plt.show()
Out [7]:
Simulation visualization [7]
Figure 3:Simulation visualization [7]

Experiment 4: Phase Space Portrait in (N1, N2)

In dynamic systems theory, a phase space portrait plots the vector field \(\left(\frac{dN_1}{dt}, \frac{dN_2}{dt}\right)\) and state trajectories starting from arbitrary initial memory conditions.

Notice that regardless of whether the cache starts empty \((0, 0)\), skewed towards light keys \((160\text{k}, 5\text{k})\), or with heavy pre-warmed keys, all trajectories converge toward the unique globally asymptotically stable equilibrium point (stable node attractor) governed by the Lyapunov stability of the linear decay terms.

python
# Generate 2D Phase Portrait
fig4, ax4 = plt.subplots(figsize=(8.5, 5.5))
apply_carbon_style(fig4, ax4)

# Grid in thousands
n1_g = np.linspace(0, 180000, 20)
n2_g = np.linspace(0, 25000, 20)
N1_m, N2_m = np.meshgrid(n1_g, n2_g)

M_m = config.s1 * N1_m + config.s2 * N2_m
ev_m = config.alpha * np.maximum(0.0, M_m - config.maxmemory)
U = 200 - N1_m / config.ttl1 - (config.s1 * N1_m / (M_m + config.eps_s)) * ev_m
V = 4 - N2_m / config.ttl2 - (config.s2 * N2_m / (M_m + config.eps_s)) * ev_m

norm_v = np.hypot(U, V)
norm_v[norm_v == 0] = 1
ax4.quiver(N1_m / 1000.0, N2_m / 1000.0, U / norm_v, V / norm_v, color='#4589ff', alpha=0.35, angles='xy', scale=25)

# Trajectories
trajectories = [
    [0, 0],
    [160000, 5000],
    [40000, 22000],
    [160000, 22000]
]
traj_colors = ['#08bdba', '#f1c21b', '#d12771', '#be95ff']
for init_s, col in zip(trajectories, traj_colors):
    sol_t = solve_ivp(cache_dynamics, [0, 8000], init_s, args=(200, 4), t_eval=np.linspace(0, 8000, 300))
    ax4.plot(sol_t.y[0] / 1000.0, sol_t.y[1] / 1000.0, color=col, lw=2, label=f'Init ({init_s[0]//1000}k, {init_s[1]//1000}k)')

# Maxmemory boundary
n1_bound = np.linspace(0, 180000, 100)
n2_bound = (config.maxmemory - config.s1 * n1_bound) / config.s2
ax4.plot(n1_bound / 1000.0, n2_bound / 1000.0, color='#da1e28', linestyle=':', lw=2, label='Memory Limit (4 GB)')

# Attractor
ax4.plot(N1_q / 1000.0, N2_q / 1000.0, 'o', color='#42be65', markersize=9, label=f'Stable Attractor ({N1_q/1000:.0f}k, {N2_q/1000:.1f}k)')

ax4.set_xlim(0, 180)
ax4.set_ylim(0, 25)
ax4.set_xlabel('Class 1 Objects: N1 (thousands)', fontsize=11)
ax4.set_ylabel('Class 2 Objects: N2 (thousands)', fontsize=11)
ax4.set_title('Phase Space Portrait (N1, N2) under Quadrupled Load', fontsize=12, pad=12)
ax4.legend(facecolor='#262626', edgecolor='#393939', labelcolor='#f4f4f4', loc='upper right', fontsize=8.5)
fig4.tight_layout()
plt.show()
Out [8]:
Simulation visualization [8]
Figure 4:Simulation visualization [8]

Conclusion: Common Sense for Dynamic Systems

The model portraits reveal how the final stationary number of heavy results changes as the load increases:

  • At \(\lambda_2 \in \{1,\dots,4\}\) → the growth is almost linear.
  • At \(\lambda_2 \approx 4,\dots,6\) → the system enters the memory saturation zone.
  • At \(\lambda_2 > 6\) → Redis enters the mode of constant evictions; the growth of \(N_2\) slows down almost to zero → thrashing mode.

The number \(N_2\) answers the question: "How many heavy results will, on average, reside in Redis at a given heavy task intensity \(\lambda_2\) once the system has stabilized?" The answer: close to 20–21 thousand.

New parameters (TTL, load, result size, memory) do not immediately appear in the metrics but after a couple of TTLs. This means after changing the configuration, you must wait several TTLs before evaluating the effect.

Yes, it has little to do with precise reality. Yes, there are a lot of caveats. But this is already a way to at least slightly tame the digital chaos. Systems live. They are not static. They flow, change, drift, get tired, behave strangely, sometimes delight, sometimes bite. And our job is not to try to force them to stand still, like photographs on the wall.

Our job is to understand that they are dynamic and to factor that into every architectural, engineering, and operational detail.

That's all for now. Beyond this — only practice, monitoring, and cautious optimism.