mtp_sim.py raw
1 #!/usr/bin/env python3
2 """
3 Limenka MTP-cadence simulation - core discrete-event model.
4
5 Models the settled single-lane 600s fork design:
6
7 e = nTime(i) - nTime(i-1), direct stamp delta (MTP-101 is used
8 for the activation gate only, not for e)
9 reward = R_full * e / 600 (uncapped, floor zero)
10 payload = S_max * e / 600 (uncapped; header+coinbase subtracted)
11 DAA = PI controller, error = 600 - e, per-block update.
12 P (proportional) = kp on the instantaneous error, symmetric.
13 I (integral) = ki on the EMA accumulator, asymmetric
14 raise-slow/lower-fast for flood resistance.
15 D (derivative) = kd on the band-passed error difference -
16 measured dead in this system, kd = 0.
17 time = fork blocks stamped <= now + F (default 60s), honest NTP
18 assumed: no stamp may exceed arrival + F
19
20 Attackers may stamp blocks up to arrival + F (reward inflation via time
21 fakery bounded by F) and self-pad blocks to the full payload allowance
22 with zero-cost self-pay txs, so realized attacker block size equals the
23 allowance, not the mempool backlog.
24
25 Scenarios: steady, ramp_slow, ramp_fast, collapse, flood, monopoly,
26 faketime, stall (see scenarios()).
27
28 Metrics: mean/median/p90 interval, issuance vs parity, difficulty
29 max/final, stall window after flood exit, reward split, backlog stats,
30 realized max block size.
31
32 Run: python3 mtp_sim.py --scenario flood --blocks 30000
33 """
34 import argparse
35 import bisect
36 import random
37
38 T = 600 # target interval, seconds
39 S_MAX = 4_000_000 # standard payload bytes (witness discount removed)
40 R_FULL = 1.0 # normalized full subsidy: 1 per 600s
41
42 PID_SCALE = 1_000_000
43 DEN_MIN = 100_000 # per-block target move clamp: 0.1x
44 DEN_MAX = 10_000_000 # per-block target move clamp: 10x
45
46
47 class World:
48 """Piecewise-linear hashrate schedule. const(t, h, a), ramp(t0, t1, h0, h1).
49
50 Segments carry the attacker's stamp mode: 'honest' (stamp at arrival),
51 'future' (stamp at arrival+F, max allowed), 'past' (stamp at MTP floor,
52 the earliest valid - can be ~window/2*T behind wall clock)."""
53
54 def __init__(self):
55 self.segs = []
56
57 def const(self, start, h, a=0.0, mode="honest"):
58 self.segs.append((start, None, h, h, a, mode))
59
60 def ramp(self, start, end, h0, h1, a=0.0, mode="honest"):
61 self.segs.append((start, end, h0, h1, a, mode))
62
63 def rates_at(self, t):
64 h = a = 0.0
65 mode = "honest"
66 for start, end, h0, h1, aa, md in self.segs:
67 if t >= start:
68 h, a, mode = h0, aa, md
69 if end is None:
70 continue
71 if t < end:
72 frac = (t - start) / (end - start)
73 h = h0 + (h1 - h0) * frac
74 else:
75 h = h1
76 return h, a, mode
77
78 def next_boundary(self, t):
79 """Next schedule change strictly after t, or None."""
80 nxt = None
81 for start, end, _, _, _, _ in self.segs:
82 for cand in (start, end):
83 if cand is None:
84 continue
85 if cand > t:
86 nxt = cand if nxt is None else min(nxt, cand)
87 return nxt
88
89
90 class Chain:
91 def __init__(self, cfg):
92 self.cfg = cfg
93 self.d = 1.0 # difficulty multiplier (1 = parent)
94 self.t = 0.0 # wall clock, seconds
95 self.mtp_sorted = [] # sorted stamps of last W blocks
96 self.mtp_deque = []
97 self.avg_error = 0.0 # EMA of error (P accumulator)
98 self.err_lp = 0.0 # low-passed error (D state)
99 self.err_lp_prev = 0.0 # previous low-passed error (D state)
100 self.last_block_t = 0.0 # completion time of last block
101 self.height = 0
102 self.issuance = 0.0
103 self.backlog = 0.0 # mempool backlog, bytes
104 self.max_backlog = 0.0
105 self.max_block_size = 0.0
106 self.honest_reward = 0.0
107 self.attacker_reward = 0.0
108 self.blocks = [] # (arrival, attacker, e, d, size)
109 self.stamp_violations = 0
110 self.orphans = 0 # blocks lost to propagation races
111 # delay floor, asymmetric gains, propagation - optional cfg attrs
112 # Naming (corrected from the swapped 2021 lineage):
113 # P (proportional) = gain on the instantaneous error
114 # I (integral) = gain on the EMA accumulator (avg_error)
115 # D (derivative) = gain on the band-passed error difference
116 self.delay = getattr(cfg, "delay", 0.0)
117 self.att_hw = getattr(cfg, "att_hw", 1.0)
118 self.prop_base = getattr(cfg, "prop_base", 0.0)
119 self.bandwidth = getattr(cfg, "bandwidth", 100e6)
120 self.honest_miners = getattr(cfg, "honest_miners", 1)
121 # DAA setpoint, separate from the reward time basis T (600s).
122 # Calibrated above 600 to cancel the asymmetric-integral
123 # rectification bias (the ~17s steady dilation).
124 self.target = getattr(cfg, "target", T)
125 self.kp = cfg.kp
126 self.ki_up = getattr(cfg, "ki_up", cfg.ki)
127 self.ki_down = getattr(cfg, "ki_down", cfg.ki)
128
129 def mtp(self):
130 n = len(self.mtp_sorted)
131 return self.mtp_sorted[(n - 1) // 2] if n else 0
132
133 def advance_time(self, t_new):
134 dt = t_new - self.t
135 if dt > 0:
136 self.backlog += self.cfg.fill * dt
137 self.max_backlog = max(self.max_backlog, self.backlog)
138 self.t = t_new
139
140 def push_stamp(self, s):
141 bisect.insort(self.mtp_sorted, s)
142 self.mtp_deque.append(s)
143 if len(self.mtp_deque) > self.cfg.window:
144 old = self.mtp_deque.pop(0)
145 del self.mtp_sorted[bisect.bisect_left(self.mtp_sorted, old)]
146
147 def add_block(self, arrival, attacker, seed=False, mode="honest"):
148 mtp_prev = self.mtp()
149 last_stamp = self.mtp_deque[-1] if self.mtp_deque else arrival
150 cap = arrival + self.cfg.future
151 if self.cfg.monotonic:
152 # strict stamp monotonicity: nTime > prev nTime (fork rule)
153 floor = last_stamp + 1
154 else:
155 # parent-style rule: nTime > MTP(prev) only
156 floor = mtp_prev + 1
157 if mode == "past":
158 # stamp at the earliest consensus-valid time, possibly far in
159 # the past relative to wall clock (MTP floor only)
160 stamp = floor
161 elif mode == "future":
162 # stamp as far ahead as the future limit allows
163 stamp = min(cap, max(floor, arrival))
164 else:
165 # honest: stamp at real arrival
166 stamp = min(cap, max(arrival, floor))
167 if cap < floor:
168 # both rules cannot be satisfied at `arrival`: the miner waits
169 # for the clock to catch up, then stamps at the floor
170 self.stamp_violations += 1
171 arrival = floor - self.cfg.future
172 cap = arrival + self.cfg.future
173 stamp = floor
174
175 if seed:
176 # first block: continue parent difficulty, no DAA update
177 self.advance_time(arrival)
178 self.last_block_t = arrival
179 self.push_stamp(stamp)
180 self.height += 1
181 return
182
183 self.advance_time(arrival)
184 self.last_block_t = arrival
185 self.push_stamp(stamp)
186 if self.cfg.ebasis == "mtp":
187 e = max(1, self.mtp() - mtp_prev)
188 else: # "direct": stamp-to-stamp delta of consecutive blocks
189 e = max(1, stamp - last_stamp)
190
191 reward = R_FULL * e / T
192 allowance = S_MAX * e / T
193 size = allowance if attacker else min(allowance, self.backlog)
194 self.backlog = max(0.0, self.backlog - size)
195 self.max_block_size = max(self.max_block_size, size)
196
197 self.issuance += reward
198 if attacker:
199 self.attacker_reward += reward
200 else:
201 self.honest_reward += reward
202
203 error = self.target - e
204 self.avg_error = self.avg_error * (1.0 - self.cfg.alpha) + error * self.cfg.alpha
205 # band-pass derivative: difference of the low-passed error.
206 # The low-pass kills high-frequency noise; the differencing
207 # removes the DC component - a band-pass in the standard PID
208 # filtered-derivative sense. (Measured dead in this system - the
209 # single-sample exponential measurement has SNR 1, the derivative
210 # carries ~3-5x more noise than signal. kd stays 0.)
211 self.err_lp_prev = self.err_lp
212 self.err_lp = self.err_lp * (1.0 - self.cfg.lp) + error * self.cfg.lp
213 deriv = self.err_lp - self.err_lp_prev
214 # P: proportional on the instantaneous error (symmetric).
215 # I: integral on the EMA accumulator, asymmetric raise-slow /
216 # lower-fast for flood resistance.
217 ki = self.ki_up if error > 0 else self.ki_down
218 correction = (ki * self.avg_error + self.kp * error +
219 self.cfg.kd * deriv) / PID_SCALE
220 denom = max(DEN_MIN, min(DEN_MAX, PID_SCALE + correction))
221 self.d *= denom / PID_SCALE
222
223 self.height += 1
224 self.blocks.append((arrival, attacker, e, self.d, size))
225
226
227 def scenarios(cfg):
228 w = World()
229 if cfg.scenario == "steady":
230 w.const(0, 1.0)
231 elif cfg.scenario == "ramp_slow":
232 w.ramp(0, 30 * 86400, 1.0, 2.0)
233 elif cfg.scenario in ("ramp_fast", "ramp100", "ramp1000", "ramp1e6"):
234 w.ramp(0, 86400, 1.0, cfg.ramp_x)
235 elif cfg.scenario in ("collapse", "collapse99"):
236 w.const(0, 1.0)
237 w.const(100 * T, 1.0 - cfg.drop)
238 elif cfg.scenario == "flood":
239 w.const(0, 1.0)
240 w.const(100 * T, 1.0, cfg.flood_x)
241 w.const(100 * T + cfg.flood_min * 60, 1.0)
242 elif cfg.scenario == "monopoly":
243 w.const(0, 0.1, 0.9)
244 elif cfg.scenario in ("faketime", "future_attack", "future_minority",
245 "future_majority"):
246 share = cfg.attack_share if cfg.scenario != "faketime" else 0.6
247 w.const(0, 1.0 - share, share, mode="future")
248 elif cfg.scenario in ("past_attack", "past_minority", "past_majority"):
249 share = cfg.attack_share
250 w.const(0, 1.0 - share, share, mode="past")
251 elif cfg.scenario == "stall":
252 w.const(0, 1.0)
253 w.const(100 * T, cfg.stall_frac)
254 w.const(100 * T + cfg.stall_days * 86400, 1.0)
255 else:
256 raise SystemExit(f"unknown scenario {cfg.scenario}")
257 return w
258
259
260 def flood_exit(cfg):
261 if cfg.scenario == "flood":
262 return 100 * T + cfg.flood_min * 60
263 return None
264
265
266 def run(cfg):
267 """Run the simulation. Returns a dict of metrics."""
268 rng = random.Random(cfg.seed)
269 world = scenarios(cfg)
270 chain = Chain(cfg)
271 exit_t = flood_exit(cfg)
272
273 chain.add_block(0.0, False, seed=True)
274 for _ in range(cfg.blocks):
275 h, a, mode = world.rates_at(chain.t)
276 total = h + a
277 if total <= 0:
278 break
279 dt = rng.expovariate(total / (T * chain.d))
280 arr = chain.t + dt
281 nxt = world.next_boundary(chain.t)
282 if nxt is not None and arr >= nxt:
283 chain.advance_time(nxt) # rate changes before any block arrives
284 continue
285 attacker = rng.random() < a / total
286 hw = chain.att_hw if attacker else 1.0
287 # sequential delay floor: block can't complete before prev + D*hw
288 min_arr = chain.last_block_t + chain.delay * hw
289 if arr < min_arr:
290 arr = min_arr
291 if nxt is not None and arr >= nxt:
292 chain.advance_time(nxt)
293 continue
294 # propagation race: competing miners finishing within our
295 # propagation window orphan our block. The delay floor serializes
296 # block production when D > propagation, so competitors cannot
297 # exist inside the window (compet_window -> 0).
298 if chain.prop_base > 0:
299 est_e = max(1.0, arr - chain.last_block_t)
300 est_size = (S_MAX * est_e / T) if attacker else min(S_MAX * est_e / T, chain.backlog)
301 prop = chain.prop_base + est_size / chain.bandwidth
302 hw_other = chain.att_hw if not attacker else 1.0
303 # competing hashrate: the rest of the honest population (split
304 # into N miners) plus the attacker side
305 if attacker:
306 other_rate = h / (T * chain.d)
307 else:
308 other_rate = (h * (chain.honest_miners - 1) / chain.honest_miners + a) / (T * chain.d)
309 # competitor completion spacing: PoW mean, floored by their
310 # delay (delay-bound competitors serialize at the floor)
311 spacing = max((1.0 / other_rate) if other_rate > 0 else float("inf"),
312 chain.delay * hw_other)
313 p_orphan = min(1.0, 2.0 * prop / spacing) * 0.5
314 if rng.random() < p_orphan:
315 chain.orphans += 1
316 continue
317 chain.add_block(arr, attacker, mode=mode if attacker else "honest")
318
319 ints = sorted(b[2] for b in chain.blocks if b[2] > 0)
320 n = max(1, len(ints))
321 parity = R_FULL * chain.t / T
322 res = {
323 "scenario": cfg.scenario,
324 "blocks": chain.height,
325 "sim_hours": chain.t / 3600,
326 "mean_int": sum(ints) / n,
327 "med_int": ints[n // 2],
328 "p90_int": ints[int(n * 0.9)],
329 "issuance_ratio": chain.issuance / parity if parity > 0 else 0.0,
330 "d_max": max(b[3] for b in chain.blocks) if chain.blocks else 1.0,
331 "d_final": chain.d,
332 "max_backlog": chain.max_backlog,
333 "max_block_size": chain.max_block_size,
334 "honest_reward": chain.honest_reward,
335 "attacker_reward": chain.attacker_reward,
336 "stamp_violations": chain.stamp_violations,
337 "stall_first_blk": float("nan"),
338 "stall_restore": float("nan"),
339 }
340
341 if exit_t is not None:
342 post = [b for b in chain.blocks if b[0] >= exit_t and not b[1]]
343 if post:
344 res["stall_first_blk"] = post[0][0] - exit_t
345 good = 0
346 for b in post:
347 if abs(b[2] - T) <= 0.5 * T:
348 good += 1
349 else:
350 good = 0
351 if good >= 3:
352 res["stall_restore"] = b[0] - exit_t
353 break
354 return res
355
356
357 def main():
358 ap = argparse.ArgumentParser(description="Limenka MTP-cadence sim")
359 ap.add_argument("--scenario", default="steady")
360 ap.add_argument("--blocks", type=int, default=30000)
361 ap.add_argument("--seed", type=int, default=42)
362 ap.add_argument("--window", type=int, default=101)
363 ap.add_argument("--future", type=int, default=60)
364 ap.add_argument("--kp", type=float, default=50e6)
365 ap.add_argument("--ki", type=float, default=10e6)
366 ap.add_argument("--kd", type=float, default=16.7e6)
367 ap.add_argument("--alpha", type=float, default=0.10)
368 ap.add_argument("--lp", type=float, default=0.15)
369 ap.add_argument("--fill", type=float, default=6600)
370 ap.add_argument("--ebasis", choices=["mtp", "direct"], default="mtp")
371 ap.add_argument("--drop", type=float, default=0.9)
372 ap.add_argument("--flood-x", type=float, default=50)
373 ap.add_argument("--flood-min", type=float, default=30)
374 ap.add_argument("--stall-frac", type=float, default=0.001)
375 ap.add_argument("--stall-days", type=float, default=4)
376 ap.add_argument("--ramp-x", type=float, default=10)
377 ap.add_argument("--attack-share", type=float, default=0.2)
378 ap.add_argument("--monotonic", action="store_true")
379 cfg = ap.parse_args()
380
381 res = run(cfg)
382 for k, v in res.items():
383 if isinstance(v, float):
384 print(f"{k:18s} {v:,.3f}")
385 else:
386 print(f"{k:18s} {v}")
387
388
389 if __name__ == "__main__":
390 main()
391