feat(scenario): scenario runner, extensions and working copy
scenario.py runs all five simulators of the core from a JSON scenario file, with sweeps, and reproduces 20 of the note's data files (14 byte for byte, 6 column-wise). Extensions live outside the core and must reproduce it exactly in a limit case; paadrag_tidsserie is the worked example. arbeidskopi.py gives the user a copy of the model they can change. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
This commit is contained in:
parent
464ad38b49
commit
03151cfb3b
SSH key fingerprint:
SHA256:JakMjO6FTBBzN0Bhfj9saOoEjaFxlSdYuZQQpM/lF9Q
36 changed files with 3580 additions and 5 deletions
8
modell/utvidelser/__init__.py
Normal file
8
modell/utvidelser/__init__.py
Normal file
|
|
@ -0,0 +1,8 @@
|
|||
# Copyright (c) 2026 Kjell Tore Guttormsen. MIT License, see LICENSE.
|
||||
"""Utvidelser: mekanismer utenfor kjernen.
|
||||
|
||||
Hver modul her har NAVN, BESKRIVELSE, PARAMETERE, simuler(scenario) og
|
||||
grensetilfelle(). Regelen er den samme som i notatet: en utvidelse skal ha et
|
||||
grensetilfelle som gjenskaper kjernen eksakt, og en test som viser det. Se
|
||||
docs/utvidelser.md.
|
||||
"""
|
||||
258
modell/utvidelser/paadrag_tidsserie.py
Normal file
258
modell/utvidelser/paadrag_tidsserie.py
Normal file
|
|
@ -0,0 +1,258 @@
|
|||
# Copyright (c) 2026 Kjell Tore Guttormsen. MIT License, see LICENSE.
|
||||
"""Utvidelse: det strukturelle pådraget følger en tidsserie.
|
||||
|
||||
I kjernen (sim.py, Omega) er pådraget Omega0 før måned t_shock og Omega0 * shock
|
||||
etterpå: ett varig sjokk. Her følger pådraget Omega0 * f(t), der f er en
|
||||
tidsserie brukeren oppgir, som trappetrinn eller stykkevis lineær. Da kan en
|
||||
forbigående bølge, et sjokk i flere trinn eller en gradvis endring kjøres.
|
||||
|
||||
Styringen er P-regelen fra kjernen. Med kjeden (parametere.CH0.Z endelig) er
|
||||
regnestegene de samme som i simulate_chain, med tiltakene i kjeden (andel
|
||||
utenom, kapasitet, trommel-buffer-tau). Uten kjeden (CH0.Z = "inf") er de de
|
||||
samme som i simulate.
|
||||
|
||||
Grensetilfelle: et trappetrinn lik notatets sjokk (faktor 1,3 fra måned 12)
|
||||
gjenskaper simulate_chain("P", 3.0, CH0) og, med CH0.Z = inf, simulate("P", 3.0)
|
||||
bit for bit. Regnestegene under er kopier av løkkene i sim.py; endres de der,
|
||||
må de endres likt her, og grensetilfellet viser om det er glemt.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import importlib
|
||||
import math
|
||||
import sys
|
||||
from pathlib import Path
|
||||
from typing import TYPE_CHECKING, Any, cast
|
||||
|
||||
import numpy as np
|
||||
|
||||
if TYPE_CHECKING:
|
||||
from types import ModuleType
|
||||
|
||||
import numpy.typing as npt
|
||||
|
||||
Arr = npt.NDArray[np.float64]
|
||||
|
||||
NAVN = "paadrag_tidsserie"
|
||||
BESKRIVELSE = ("det strukturelle pådraget Omega følger en tidsserie (trappetrinn eller stykkevis "
|
||||
"lineær) i stedet for ett varig sjokk; P-regel, med eller uten straffesakskjeden")
|
||||
PARAMETERE: dict[str, object] = {
|
||||
"form": "trapp", # "trapp" eller "lineaer"
|
||||
"punkter": [[12.0, 1.3]], # [måned, faktor på Omega0]; før første punkt er faktoren 1
|
||||
}
|
||||
BRUKER = ("P", "CH0") # kjernens dicts som kan overstyres under parametere
|
||||
FORMER = ("trapp", "lineaer")
|
||||
MAKS_PUNKTER = 1000
|
||||
KJERNE = Path(__file__).resolve().parents[1] / "kjerne"
|
||||
|
||||
|
||||
def _sim() -> ModuleType:
|
||||
"""Kjernen (sim.py) fra modellmappa denne utvidelsen ligger i."""
|
||||
sim = sys.modules.get("sim")
|
||||
if sim is not None and Path(str(getattr(sim, "__file__", ""))).resolve().parent == KJERNE:
|
||||
return sim
|
||||
sys.modules.pop("sim", None)
|
||||
if str(KJERNE) in sys.path:
|
||||
sys.path.remove(str(KJERNE))
|
||||
sys.path.insert(0, str(KJERNE))
|
||||
return importlib.import_module("sim")
|
||||
|
||||
|
||||
def _punkter(raa: object) -> list[tuple[float, float]]:
|
||||
"""Valider punktene: [[måned, faktor], ...] med stigende måneder."""
|
||||
hvor = "utvidelse.punkter"
|
||||
if not isinstance(raa, list) or not 1 <= len(raa) <= MAKS_PUNKTER:
|
||||
melding = f"{hvor}: må være en liste med 1 til {MAKS_PUNKTER} punkter [måned, faktor]."
|
||||
raise ValueError(melding)
|
||||
ut: list[tuple[float, float]] = []
|
||||
for nr, p in enumerate(raa):
|
||||
ok = (isinstance(p, list) and len(p) == 2 # noqa: PLR2004 - et punkt er et par
|
||||
and all(isinstance(x, (int, float)) and not isinstance(x, bool) for x in p))
|
||||
if not ok:
|
||||
melding = f"{hvor}[{nr}]: må være et par [måned, faktor] med tall, ikke {p!r}."
|
||||
raise ValueError(melding)
|
||||
tid, faktor = float(p[0]), float(p[1])
|
||||
if not (math.isfinite(tid) and math.isfinite(faktor)) or tid < 0 or faktor < 0:
|
||||
melding = f"{hvor}[{nr}]: måned og faktor må være tall fra og med 0."
|
||||
raise ValueError(melding)
|
||||
if ut and tid <= ut[-1][0]:
|
||||
melding = f"{hvor}[{nr}]: månedene må stige (måned {tid:g} kommer etter {ut[-1][0]:g})."
|
||||
raise ValueError(melding)
|
||||
ut.append((tid, faktor))
|
||||
return ut
|
||||
|
||||
|
||||
def faktor(tid: float, form: str, punkter: list[tuple[float, float]]) -> float:
|
||||
"""f(t): 1 før første punkt; trapp holder verdien til neste punkt, lineaer interpolerer."""
|
||||
if tid < punkter[0][0]:
|
||||
return 1.0
|
||||
if form == "trapp":
|
||||
verdi = 1.0
|
||||
for t0, f0 in punkter:
|
||||
if tid >= t0:
|
||||
verdi = f0
|
||||
else:
|
||||
break
|
||||
return verdi
|
||||
for (t0, f0), (t1, f1) in zip(punkter, punkter[1:]):
|
||||
if t0 <= tid < t1:
|
||||
return f0 + (f1 - f0) * (tid - t0) / (t1 - t0)
|
||||
return punkter[-1][1]
|
||||
|
||||
|
||||
def sjekk(scenario: dict[str, object]) -> None:
|
||||
"""Kalles av scenariokjøreren før kjøringen; ValueError gir en norsk feilmelding."""
|
||||
sim = _sim()
|
||||
styring = cast("dict[str, Any]", scenario["styring"])
|
||||
if styring["regulator"] != "P":
|
||||
melding = 'styring.regulator: utvidelsen paadrag_tidsserie bruker bare "P".'
|
||||
raise ValueError(melding)
|
||||
par = cast("dict[str, dict[str, Any]]", scenario["parametere"])
|
||||
for nokkel in ("shock", "t_shock"):
|
||||
if par["P"][nokkel] != sim.P[nokkel]:
|
||||
melding = (f"parametere.P.{nokkel}: brukes ikke av paadrag_tidsserie. Pådraget følger "
|
||||
"utvidelse.punkter.")
|
||||
raise ValueError(melding)
|
||||
if math.isinf(par["CH0"]["Z"]) and scenario.get("tiltak"):
|
||||
melding = "tiltak: virker bare i kjeden, og parametere.CH0.Z er inf (ingen kjede)."
|
||||
raise ValueError(melding)
|
||||
u = cast("dict[str, Any]", scenario["utvidelse"])
|
||||
if u["form"] not in FORMER:
|
||||
melding = f'utvidelse.form: må være "trapp" eller "lineaer", ikke {u["form"]!r}.'
|
||||
raise ValueError(melding)
|
||||
_punkter(u["punkter"])
|
||||
|
||||
|
||||
def _deriv(tilstand: tuple[float, float, float, float], om: float,
|
||||
p: dict[str, Any]) -> tuple[float, float, float]:
|
||||
"""Som deriv i sim.py, med pådraget om i stedet for Omega(t, p). Tilstand: K, U, L, T."""
|
||||
sim = _sim()
|
||||
krim, trygg, tillit, innsats = tilstand
|
||||
d_krim = om - sim.removal(krim, tillit, innsats, p) - p["delta"] * krim
|
||||
u_eq = np.clip(p["U0"] - p["psiU"] * (krim - 1.0), 0.0, 1.0)
|
||||
d_trygg = p["kU"] * (u_eq - trygg)
|
||||
over = max(0.0, innsats - p["Tprop"] * krim)
|
||||
d_tillit = (p["kL"] * (p["Lam0"] + p["mu1"] * (trygg - p["U0"]) - tillit)
|
||||
- p["mu2"] * tillit * over ** 2)
|
||||
return d_krim, d_trygg, d_tillit
|
||||
|
||||
|
||||
def _uten_kjede(t: Arr, om: Arr, styring: tuple[float, float], p: dict[str, Any]) -> dict[str, object]:
|
||||
"""Løkka i simulate (P-regelen), med pådraget fra tidsserien. Styring: (tau, kp)."""
|
||||
tau, kp = styring
|
||||
sim = _sim()
|
||||
dt = sim.DT
|
||||
n = len(t)
|
||||
krim = np.ones(n)
|
||||
trygg = np.full(n, p["U0"])
|
||||
tillit = np.full(n, p["Lam0"])
|
||||
innsats = np.full(n, p["T0"])
|
||||
lag = round(tau / dt)
|
||||
for k in range(n - 1):
|
||||
k_obs = krim[k - lag] if k - lag >= 0 else 1.0
|
||||
innsats[k] = np.clip(p["T0"] + kp * (k_obs - p["Kstar"]), 0.0, p["Tmax"])
|
||||
d_krim, d_trygg, d_tillit = _deriv((krim[k], trygg[k], tillit[k], innsats[k]), om[k], p)
|
||||
krim[k + 1] = krim[k] + dt * d_krim
|
||||
trygg[k + 1] = trygg[k] + dt * d_trygg
|
||||
tillit[k + 1] = np.clip(tillit[k] + dt * d_tillit, 1e-4, 1 - 1e-4)
|
||||
innsats[-1] = innsats[-2]
|
||||
return {"t": t, "K": krim, "U": trygg, "L": tillit, "T": innsats}
|
||||
|
||||
|
||||
def _med_kjede(t: Arr, om: Arr, styring: tuple[float, float], p: dict[str, Any],
|
||||
ch: dict[str, Any]) -> dict[str, object]:
|
||||
"""Løkka i simulate_chain (P-regelen og trommel-buffer-tau), med pådraget fra tidsserien."""
|
||||
tau, kp = styring
|
||||
sim = _sim()
|
||||
dt = sim.DT
|
||||
n = len(t)
|
||||
krim = np.ones(n)
|
||||
trygg = np.full(n, p["U0"])
|
||||
tillit = np.full(n, p["Lam0"])
|
||||
innsats = np.full(n, p["T0"])
|
||||
ko, tid_d, inn_d, ut_d = np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n)
|
||||
q0, tau_ref = sim.init_chain(ch, p)
|
||||
ko[0] = q0
|
||||
lag = round(tau / dt)
|
||||
for k in range(n - 1):
|
||||
k_obs = krim[k - lag] if k - lag >= 0 else 1.0
|
||||
innsats[k] = np.clip(p["T0"] + kp * (k_obs - p["Kstar"]), 0.0, p["Tmax"])
|
||||
if ch["dbr"] is not None:
|
||||
z, xi = sim.chain_at(ch, t[k])
|
||||
t_tak = ch["dbr"] * z / (xi * p["gamma"] * tillit[k] ** p["eta"] * k_obs)
|
||||
innsats[k] = min(innsats[k], max(t_tak, 0.0))
|
||||
d_krim, d_trygg, d_tillit, d_ko, tid_d[k], inn_d[k], ut_d[k] = sim.deriv_chain(
|
||||
krim[k], trygg[k], tillit[k], innsats[k], ko[k], t[k], p, ch, tau_ref, Om0=om[k])
|
||||
krim[k + 1] = krim[k] + dt * d_krim
|
||||
trygg[k + 1] = trygg[k] + dt * d_trygg
|
||||
tillit[k + 1] = np.clip(tillit[k] + dt * d_tillit, 1e-4, 1 - 1e-4)
|
||||
ko[k + 1] = max(ko[k] + dt * d_ko, 0.0)
|
||||
innsats[-1] = innsats[-2]
|
||||
tid_d[-1] = tid_d[-2]
|
||||
inn_d[-1] = inn_d[-2]
|
||||
ut_d[-1] = ut_d[-2]
|
||||
par = [sim.chain_at(ch, float(x)) for x in t]
|
||||
z_t = np.array([float(a) for a, _ in par])
|
||||
return {"t": t, "K": krim, "U": trygg, "L": tillit, "T": innsats, "Q": ko, "tauD": tid_d,
|
||||
"JD": inn_d, "RD": ut_d, "Z": z_t, "Z_rel": z_t / ch["Z"],
|
||||
"xi": np.array([float(b) for _, b in par]), "tau_ref": tau_ref}
|
||||
|
||||
|
||||
def simuler(scenario: dict[str, object]) -> dict[str, object]:
|
||||
"""Kjør modellen med pådraget fra tidsserien. Nøkler som i kjernen (t, K, U, L, T, ...)."""
|
||||
sim = _sim()
|
||||
styring = cast("dict[str, Any]", scenario["styring"])
|
||||
par = cast("dict[str, dict[str, Any]]", scenario["parametere"])
|
||||
u = cast("dict[str, Any]", scenario["utvidelse"])
|
||||
p, ch = par["P"], par["CH0"]
|
||||
punkter = _punkter(u["punkter"])
|
||||
slutt = float(cast("float", scenario["slutt"]))
|
||||
n = int(slutt / sim.DT) + 1
|
||||
t = np.arange(n) * sim.DT
|
||||
f = np.array([faktor(x, u["form"], punkter) for x in t])
|
||||
om = np.array([p["Omega0"] * x for x in f])
|
||||
tau_kp = (float(styring["forsinkelse"]), float(styring["reaktivitet"]))
|
||||
ut = _uten_kjede(t, om, tau_kp, p) if math.isinf(ch["Z"]) else _med_kjede(t, om, tau_kp, p, ch)
|
||||
ut["Omega_faktor"] = f
|
||||
return ut
|
||||
|
||||
|
||||
def _scenario(p: dict[str, Any], ch: dict[str, Any], slutt: float) -> dict[str, object]:
|
||||
return {
|
||||
"styring": {"regulator": "P", "forsinkelse": 3.0, "reaktivitet": 3.0},
|
||||
"slutt": slutt,
|
||||
"parametere": {"P": p, "CH0": ch},
|
||||
"tiltak": {},
|
||||
"utvidelse": {"form": "trapp", "punkter": [[float(p["t_shock"]), float(p["shock"])]]},
|
||||
}
|
||||
|
||||
|
||||
def _avvik(egen: dict[str, object], fasit: dict[str, Any], navn: tuple[str, ...]) -> float:
|
||||
"""Største absolutte avvik over tidsseriene (inf hvis lengdene er ulike)."""
|
||||
storst = 0.0
|
||||
for n in navn:
|
||||
a, b = np.asarray(egen[n]), np.asarray(fasit[n])
|
||||
if a.shape != b.shape:
|
||||
return math.inf
|
||||
storst = max(storst, float(np.abs(a - b).max()))
|
||||
return storst
|
||||
|
||||
|
||||
def grensetilfelle() -> list[tuple[str, float]]:
|
||||
"""Trappetrinn lik notatets sjokk gjenskaper kjernen; avvikene skal være eksakt 0."""
|
||||
sim = _sim()
|
||||
slutt = float(sim.TEND)
|
||||
egen = simuler(_scenario(dict(sim.P), dict(sim.CH0), slutt))
|
||||
fasit = sim.simulate_chain("P", 3.0, sim.CH0)
|
||||
d_kjede = max(_avvik(egen, fasit, ("t", "K", "U", "L", "T", "Q", "tauD", "JD", "RD")),
|
||||
abs(float(cast("float", egen["tau_ref"])) - float(fasit["tau_ref"])))
|
||||
egen0 = simuler(_scenario(dict(sim.P), dict(sim.CH0, Z=math.inf), slutt))
|
||||
t, krim, trygg, tillit, innsats = sim.simulate("P", 3.0)
|
||||
d_uten = _avvik(egen0, {"t": t, "K": krim, "U": trygg, "L": tillit, "T": innsats},
|
||||
("t", "K", "U", "L", "T"))
|
||||
return [
|
||||
('paadrag_tidsserie: trapp lik notatets sjokk, med kjeden = simulate_chain("P", 3, CH0)',
|
||||
d_kjede),
|
||||
('paadrag_tidsserie: trapp lik notatets sjokk, CH0.Z = inf = simulate("P", 3)', d_uten),
|
||||
]
|
||||
Loading…
Add table
Add a link
Reference in a new issue