Kjell Tore Guttormsen
464ad38b49
The core in modell/kjerne/ is an unmodified copy of sim.py, diagnose.py and forlop.py from den-dynamiske-rettsstaten @2badc25, pinned by sha256. kjor.py is the single entry point; kontroll.py reproduces the note's 22 data files and printed results. Skills: diagnose and forlop (Norwegian). Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
1077 lines
57 KiB
Python
1077 lines
57 KiB
Python
"""Simulering av den aggregerte modellen i artikkelen 'Den dynamiske rettsstaten'
|
|
(Guttormsen 2026), avsnitt 8.
|
|
|
|
Tidsenhet: måned. Alle nivåer normalisert slik at K = 1 er referansenivået.
|
|
Skriver CSV-filer som pgfplots leser direkte, og skriver ut nøkkeltallene som
|
|
brukes i teksten (tabell over parametere, avsnitt 8, regneeksemplene i
|
|
avsnitt 7.6 og 9.3).
|
|
|
|
Del 1: styringsregimene A-E uten straffesakskjede (fig:sim, fig:stab).
|
|
Del 2: den koblede modellen med domstolskø (avsnitt 7.3 og 8.2, fig:kjedesim),
|
|
inkludert grensetilfelle-kontroll, scenarioer, følsomhet og MPC.
|
|
Del 3: robusthet (avsnitt 8.3 og 8.4): politiets kapasitet (Holling type II,
|
|
lign. for F_sat), prioritert tokjassekø og flerserver-effekt i domstolen,
|
|
og selvforsterkende rekruttering (vippepunkt). Alle utvidelser har
|
|
grensetilfeller som gjenskaper Del 1 og 2 eksakt; kontrolltallene skrives ut.
|
|
Del 4: endogene tiltak (avsnitt 8.5): underordning xi_D, kapasitet Z_D og
|
|
prosessuell kvalitet q bestemmes av beslutningsregler som svarer på
|
|
observasjoner med forsinkelse (faste vedtak, desentraliserte regler og en
|
|
koordinert regel etter metoden i avsnitt 7.7). Faste tidsplaner gjenskaper
|
|
scenario 2-5 eksakt.
|
|
"""
|
|
import math
|
|
import numpy as np
|
|
from scipy.optimize import brentq, minimize_scalar
|
|
|
|
# ---------------- parametere (stiliserte, se tabell i artikkelen) -------------
|
|
P = dict(
|
|
gamma=0.125, # effektivitet per enhet tiltak (gamma * Lambda^eta * T = 0.1 ved likevekt)
|
|
eta=1.0, # tillitselastisitet
|
|
delta=0.02, # naturlig frafall per måned
|
|
Omega0=0.12, # strukturelt pådrag ved referanse (gir K=1 ved T=1, Lambda=0.8)
|
|
shock=1.30, # makrosjokk: Omega øker 30 % ...
|
|
t_shock=12.0, # ... fra måned 12
|
|
T0=1.0, Tmax=3.0, Kstar=1.0,
|
|
# trygghet og tillit
|
|
kU=0.25, U0=0.80, psiU=0.30,
|
|
kL=0.10, mu1=0.50, mu2=0.30, Tprop=1.6, # T_prop(K) = Tprop * K
|
|
Lam0=0.80,
|
|
ksat=np.inf, # politiets halvmetningsnivå K_sat (Holling type II); inf = ingen metning
|
|
)
|
|
DT = 0.05
|
|
TEND = 120.0
|
|
|
|
|
|
def Omega(t, p):
|
|
return p["Omega0"] * (p["shock"] if t >= p["t_shock"] else 1.0)
|
|
|
|
|
|
def removal(K, L, T, p):
|
|
"""Fjerningsleddet F = gamma' Lambda^eta T K / (1 + K/K_sat), der
|
|
gamma' = gamma (1 + K_ref/K_sat) er rekalibrert slik at F er uendret i
|
|
referanselikevekten (K_ref = 1). K_sat = inf gir eksakt lign. (12)."""
|
|
g_eff = p["gamma"] * (1.0 + 1.0 / p["ksat"])
|
|
return g_eff * L ** p["eta"] * T * K / (1.0 + K / p["ksat"])
|
|
|
|
|
|
def deriv(K, U, L, T, t, p):
|
|
dK = Omega(t, p) - removal(K, L, T, p) - p["delta"] * K
|
|
Ueq = np.clip(p["U0"] - p["psiU"] * (K - 1.0), 0.0, 1.0)
|
|
dU = p["kU"] * (Ueq - U)
|
|
over = max(0.0, T - p["Tprop"] * K)
|
|
dL = p["kL"] * (p["Lam0"] + p["mu1"] * (U - p["U0"]) - L) - p["mu2"] * L * over ** 2
|
|
return dK, dU, dL
|
|
|
|
|
|
def simulate(controller, tau, p=P, kp=3.0):
|
|
n = int(TEND / DT) + 1
|
|
t = np.arange(n) * DT
|
|
K = np.ones(n); U = np.full(n, p["U0"]); L = np.full(n, p["Lam0"]); T = np.full(n, p["T0"])
|
|
lag = int(round(tau / DT))
|
|
Tplan = p["T0"]
|
|
for k in range(n - 1):
|
|
Kobs = K[k - lag] if k - lag >= 0 else 1.0 # forsinket observasjon
|
|
if controller == "P":
|
|
T[k] = np.clip(p["T0"] + kp * (Kobs - p["Kstar"]), 0.0, p["Tmax"])
|
|
elif controller == "MPC":
|
|
# beslutning tas én gang per måned (budsjett-/planleggingssyklus)
|
|
if abs((t[k] % 1.0)) < 1e-9 or k == 0:
|
|
Tplan = mpc_step(k, K, U, L, T, t, lag, p)
|
|
T[k] = Tplan
|
|
dK, dU, dL = deriv(K[k], U[k], L[k], T[k], t[k], p)
|
|
K[k + 1] = K[k] + DT * dK
|
|
U[k + 1] = U[k] + DT * dU
|
|
L[k + 1] = np.clip(L[k] + DT * dL, 1e-4, 1 - 1e-4)
|
|
T[-1] = T[-2]
|
|
return t, K, U, L, T
|
|
|
|
|
|
def mpc_step(k, K, U, L, T, t, lag, p, H=24.0, lamT=0.05, lamL=20.0):
|
|
"""Forenklet MPC: (1) prediker dagens tilstand fra forsinket observasjon og
|
|
kjente tidligere tiltak, (2) velg ett konstant tiltaksnivå over horisonten H
|
|
som minimerer kvadratisk tap, (3) iverksett, (4) gjenta neste måned.
|
|
Myndighetene kjenner ikke sjokket før det er observert (pådraget antas lik
|
|
siste observerte verdi)."""
|
|
i0 = max(k - lag, 0)
|
|
Kh, Uh, Lh = K[i0], U[i0], L[i0]
|
|
# anslå pådraget fra det som er observert (ingen fremtidskunnskap om sjokket)
|
|
t_obs = t[i0]
|
|
def Om(_t):
|
|
return Omega(t_obs, p)
|
|
# (1) prediksjon fram til nå
|
|
for j in range(i0, k):
|
|
dK = Om(t[j]) - p["gamma"] * Lh ** p["eta"] * T[j] * Kh - p["delta"] * Kh
|
|
Ueq = np.clip(p["U0"] - p["psiU"] * (Kh - 1.0), 0, 1)
|
|
dU = p["kU"] * (Ueq - Uh)
|
|
over = max(0.0, T[j] - p["Tprop"] * Kh)
|
|
dL = p["kL"] * (p["Lam0"] + p["mu1"] * (Uh - p["U0"]) - Lh) - p["mu2"] * Lh * over ** 2
|
|
Kh += DT * dK; Uh += DT * dU; Lh = np.clip(Lh + DT * dL, 1e-4, 1 - 1e-4)
|
|
|
|
def cost(Tc):
|
|
Kc, Uc, Lc, J = Kh, Uh, Lh, 0.0
|
|
for j in range(int(H / DT)):
|
|
dK = Om(0) - p["gamma"] * Lc ** p["eta"] * Tc * Kc - p["delta"] * Kc
|
|
Ueq = np.clip(p["U0"] - p["psiU"] * (Kc - 1.0), 0, 1)
|
|
dU = p["kU"] * (Ueq - Uc)
|
|
over = max(0.0, Tc - p["Tprop"] * Kc)
|
|
dL = p["kL"] * (p["Lam0"] + p["mu1"] * (Uc - p["U0"]) - Lc) - p["mu2"] * Lc * over ** 2
|
|
Kc += DT * dK; Uc += DT * dU; Lc = np.clip(Lc + DT * dL, 1e-4, 1 - 1e-4)
|
|
J += DT * ((Kc - p["Kstar"]) ** 2 + lamT * (Tc - p["T0"]) ** 2
|
|
+ lamL * max(0.0, p["Lam0"] - 0.05 - Lc) ** 2)
|
|
return J
|
|
|
|
res = minimize_scalar(cost, bounds=(0.0, p["Tmax"]), method="bounded",
|
|
options={"xatol": 1e-3})
|
|
return float(res.x)
|
|
|
|
|
|
def critical_delay(a, b):
|
|
"""Hayes (1950): x' = -a x - b x(t-tau). Stabil for alle tau hvis b <= a."""
|
|
if b <= a:
|
|
return np.inf
|
|
return np.arccos(-a / b) / np.sqrt(b ** 2 - a ** 2)
|
|
|
|
|
|
def write_csv(name, cols, header, every=10):
|
|
arr = np.column_stack(cols)[::every]
|
|
np.savetxt(name, arr, delimiter=" ", header=header, comments="", fmt="%.5f")
|
|
|
|
|
|
# =============================================================================
|
|
# Del 2: straffesakskjeden koblet til dynamikken (avsnitt 7.3 og 8.2)
|
|
# =============================================================================
|
|
# Domstolskø Q_D (normaliserte enheter), avgjorte saker R_D = Z_D Q_D/(Q_D + m_D Z_D),
|
|
# saksbehandlingstid tau_D = m_D + Q_D/Z_D, frafall varsigma_D Q_D (henleggelse,
|
|
# foreldelse). Tilbakekoblinger: (a) celerity: Omega_eff = Omega (1 + chi_D g),
|
|
# (b) inkapasitering: andelen vartheta av fjerningseffekten for domstolssaker
|
|
# realiseres først ved avgjørelse, (c) tillit: -mu4 g i likevektsnivået for Lambda,
|
|
# med g = D/(D + tau_h), D = max(0, tau_D - tau_D^ref).
|
|
UPS0 = 0.95 # utnyttelse før sjokket (som i avsnitt 7.6)
|
|
TAU_REF = 3.0 # saksbehandlingstid i utgangslikevekten (mnd)
|
|
ZD0 = P["gamma"] * P["Lam0"] * P["T0"] * 1.0 / UPS0 # = 0.10/0.95
|
|
VARSIGMA = 0.01 # frafall fra køen per måned
|
|
# m_D kalibreres slik at tau_D = TAU_REF i utgangslikevekten (lukket form, se tekst)
|
|
J0 = P["gamma"] * P["Lam0"] * P["T0"] * 1.0
|
|
MD = TAU_REF - J0 / (ZD0 * (1.0 / TAU_REF + VARSIGMA))
|
|
|
|
CH0 = dict(Z=ZD0, m=MD, theta=0.5, chi=0.15, mu4=0.12, tauh=9.0, varsigma=VARSIGMA,
|
|
xi=1.0, Zfun=None, xifun=None, dbr=None)
|
|
|
|
|
|
def chain_at(ch, t):
|
|
Z = ch["Zfun"](t) if ch["Zfun"] else ch["Z"]
|
|
xi = ch["xifun"](t) if ch["xifun"] else ch["xi"]
|
|
return Z, xi
|
|
|
|
|
|
def court(ch, Q, JD, t, tau_ref):
|
|
"""Returnerer (R_D, tau_D, dQ/dt, g)."""
|
|
Z, _ = chain_at(ch, t)
|
|
if Z is None or np.isinf(Z): # ingen flaskehals: R_D = J_D, lign. (12) eksakt
|
|
RD, tauD, dQ = JD, tau_ref, 0.0
|
|
else:
|
|
RD = Z * Q / (Q + ch["m"] * Z)
|
|
tauD = ch["m"] + Q / Z
|
|
dQ = JD - RD - ch["varsigma"] * Q
|
|
D = max(0.0, tauD - tau_ref)
|
|
g = D / (D + ch["tauh"]) if D > 0 else 0.0
|
|
return RD, tauD, dQ, g
|
|
|
|
|
|
def deriv_chain(K, U, L, T, Q, t, p, ch, tau_ref, Om0=None):
|
|
F = removal(K, L, T, p)
|
|
_, xi = chain_at(ch, t)
|
|
JD = xi * F
|
|
RD, tauD, dQ, g = court(ch, Q, JD, t, tau_ref)
|
|
Om = (Omega(t, p) if Om0 is None else Om0) * (1.0 + ch["chi"] * g)
|
|
dK = Om - (1.0 - xi * ch["theta"]) * F - ch["theta"] * RD - p["delta"] * K
|
|
Ueq = np.clip(p["U0"] - p["psiU"] * (K - 1.0), 0.0, 1.0)
|
|
dU = p["kU"] * (Ueq - U)
|
|
over = max(0.0, T - p["Tprop"] * K)
|
|
dL = p["kL"] * (p["Lam0"] + p["mu1"] * (U - p["U0"]) - ch["mu4"] * g - L) - p["mu2"] * L * over ** 2
|
|
return dK, dU, dL, dQ, tauD, JD, RD
|
|
|
|
|
|
def init_chain(ch, p):
|
|
"""Stasjonær kø i utgangslikevekten (K=T=1, Lambda=Lam0) og referansetiden."""
|
|
Z, xi = chain_at(ch, 0.0)
|
|
JD0 = xi * p["gamma"] * p["Lam0"] * p["T0"] * 1.0
|
|
if Z is None or np.isinf(Z):
|
|
return 0.0, 0.0
|
|
f = lambda Q: Z * Q / (Q + ch["m"] * Z) + ch["varsigma"] * Q - JD0
|
|
Q0 = brentq(f, 0.0, 1e4)
|
|
return Q0, ch["m"] + Q0 / Z
|
|
|
|
|
|
def simulate_chain(controller, tau, ch, p=P, kp=3.0, mpc_knows_chain=True):
|
|
n = int(TEND / DT) + 1
|
|
t = np.arange(n) * DT
|
|
K = np.ones(n); U = np.full(n, p["U0"]); L = np.full(n, p["Lam0"]); T = np.full(n, p["T0"])
|
|
Q = np.zeros(n); tauD = np.zeros(n); JD = np.zeros(n); RD = np.zeros(n)
|
|
Q0, tau_ref = init_chain(ch, p)
|
|
Q[0] = Q0
|
|
lag = int(round(tau / DT))
|
|
Tplan = p["T0"]
|
|
for k in range(n - 1):
|
|
Kobs = K[k - lag] if k - lag >= 0 else 1.0
|
|
if controller == "P":
|
|
T[k] = np.clip(p["T0"] + kp * (Kobs - p["Kstar"]), 0.0, p["Tmax"])
|
|
elif controller == "MPC":
|
|
if abs((t[k] % 1.0)) < 1e-9 or k == 0:
|
|
if mpc_knows_chain:
|
|
Tplan = mpc_step_chain(k, K, U, L, T, Q, t, lag, p, ch, tau_ref)
|
|
else:
|
|
Tplan = mpc_step(k, K, U, L, T, t, lag, p)
|
|
T[k] = Tplan
|
|
if ch["dbr"] is not None:
|
|
# trommel-buffer-tau: begrens innsatsen slik at J_D <= ups_max Z_D
|
|
Z, xi = chain_at(ch, t[k])
|
|
Tcap = ch["dbr"] * Z / (xi * p["gamma"] * L[k] ** p["eta"] * Kobs)
|
|
T[k] = min(T[k], max(Tcap, 0.0))
|
|
dK, dU, dL, dQ, tauD[k], JD[k], RD[k] = deriv_chain(K[k], U[k], L[k], T[k], Q[k], t[k], p, ch, tau_ref)
|
|
K[k + 1] = K[k] + DT * dK
|
|
U[k + 1] = U[k] + DT * dU
|
|
L[k + 1] = np.clip(L[k] + DT * dL, 1e-4, 1 - 1e-4)
|
|
Q[k + 1] = max(Q[k] + DT * dQ, 0.0)
|
|
T[-1] = T[-2]; tauD[-1] = tauD[-2]; JD[-1] = JD[-2]; RD[-1] = RD[-2]
|
|
return dict(t=t, K=K, U=U, L=L, T=T, Q=Q, tauD=tauD, JD=JD, RD=RD, tau_ref=tau_ref)
|
|
|
|
|
|
def mpc_step_chain(k, K, U, L, T, Q, t, lag, p, ch, tau_ref, H=24.0, lamT=0.05, lamL=20.0):
|
|
"""Samme forenklede MPC som mpc_step, men med den koblede modellen (kø,
|
|
forsinket inkapasitering, celerity og tillitsvirkning) som intern modell."""
|
|
i0 = max(k - lag, 0)
|
|
Kh, Uh, Lh, Qh = K[i0], U[i0], L[i0], Q[i0]
|
|
Om_obs = Omega(t[i0], p)
|
|
|
|
def step(Kc, Uc, Lc, Qc, Tc, tt):
|
|
dK, dU, dL, dQ, _, _, _ = deriv_chain(Kc, Uc, Lc, Tc, Qc, tt, p, ch, tau_ref, Om0=Om_obs)
|
|
return (Kc + DT * dK, Uc + DT * dU, np.clip(Lc + DT * dL, 1e-4, 1 - 1e-4),
|
|
max(Qc + DT * dQ, 0.0))
|
|
|
|
for j in range(i0, k):
|
|
Kh, Uh, Lh, Qh = step(Kh, Uh, Lh, Qh, T[j], t[j])
|
|
|
|
def cost(Tc):
|
|
Kc, Uc, Lc, Qc, J = Kh, Uh, Lh, Qh, 0.0
|
|
for j in range(int(H / DT)):
|
|
Kc, Uc, Lc, Qc = step(Kc, Uc, Lc, Qc, Tc, t[k] + j * DT)
|
|
J += DT * ((Kc - p["Kstar"]) ** 2 + lamT * (Tc - p["T0"]) ** 2
|
|
+ lamL * max(0.0, p["Lam0"] - 0.05 - Lc) ** 2)
|
|
return J
|
|
|
|
res = minimize_scalar(cost, bounds=(0.0, p["Tmax"]), method="bounded",
|
|
options={"xatol": 1e-3})
|
|
return float(res.x)
|
|
|
|
|
|
def summary(name, r):
|
|
t, K, L, T, tauD = r["t"], r["K"], r["L"], r["T"], r["tauD"]
|
|
late = t > 60
|
|
ups_end = r["JD"][-1] / r["RD"][-1] if r["RD"][-1] > 0 else float("nan")
|
|
print(f" {name:34s} Kmax={K.max():.3f} Kend={K[-1]:.3f} Tmax={T.max():.2f} Tend={T[-1]:.2f} "
|
|
f"Lmin={L.min():.3f} Lend={L[-1]:.3f} tauD(end)={tauD[-1]:5.1f} tauDmax={tauD.max():5.1f} "
|
|
f"Qend={r['Q'][-1]:.2f} ups_end={ups_end:.2f} amp(K,t>60)={K[late].max() - K[late].min():.3f}")
|
|
|
|
|
|
def extrema(r, every=20):
|
|
from scipy.signal import argrelextrema
|
|
Kk = r["K"][::every]; tk = r["t"][::every]
|
|
mx = argrelextrema(Kk, np.greater)[0]; mn = argrelextrema(Kk, np.less)[0]
|
|
return [(int(round(tk[i])), round(float(Kk[i]), 3)) for i in mx if tk[i] > 15], \
|
|
[(int(round(tk[i])), round(float(Kk[i]), 3)) for i in mn if tk[i] > 15]
|
|
|
|
|
|
# =============================================================================
|
|
# Del 3: robusthet (avsnitt 8.3 og 8.4)
|
|
# =============================================================================
|
|
# (a) Politiets kapasitet: removal() ovenfor med endelig K_sat.
|
|
# (b) Domstol med to sakstyper (A: alvorlige, B: øvrige), prioritet og c "spor".
|
|
# Kapasitetsfunksjonen for c parallelle spor bygger på M/M/c-sammenhengen
|
|
# mellom beholdning L (antall saker) og utnyttelse rho (Erlang C), invertert
|
|
# numerisk; c = 1 gir eksakt R = Z Q/(Q + m Z) som i lign. (28).
|
|
# (c) Rekruttering/smitte: Omega_soc = eps K^p / (K_h^p + K^p) i tilstrømningen.
|
|
|
|
def erlang_L(rho, c):
|
|
"""Forventet antall i systemet for M/M/c ved utnyttelse rho (Erlang C)."""
|
|
A = c * rho
|
|
terms = np.array([A ** k / math.factorial(k) for k in range(c)])
|
|
last = A ** c / math.factorial(c) / (1.0 - rho)
|
|
C = last / (terms.sum() + last)
|
|
return A + C * rho / (1.0 - rho)
|
|
|
|
|
|
_ERL = {}
|
|
|
|
|
|
def rho_of_L(L, c):
|
|
"""Inverterer L(rho) for M/M/c; c=1 lukket form rho = L/(1+L)."""
|
|
if c == 1:
|
|
return L / (1.0 + L)
|
|
if c not in _ERL:
|
|
rhos = 1.0 - np.logspace(0, -7, 1400) # fra 0 til 1 - 1e-7
|
|
rhos = np.concatenate(([0.0], rhos))
|
|
Ls = np.array([0.0] + [erlang_L(r, c) for r in rhos[1:]])
|
|
_ERL[c] = (Ls, rhos)
|
|
Ls, rhos = _ERL[c]
|
|
return float(np.interp(L, Ls, rhos, right=1.0))
|
|
|
|
|
|
def capfun(Z, Q, m, c):
|
|
"""Avgjorte saker per måned med kapasitet Z, kø Q, behandlingstid m, c spor."""
|
|
if Z <= 0.0 or Q <= 0.0:
|
|
return 0.0
|
|
if c == 1:
|
|
return Z * Q / (Q + m * Z)
|
|
return Z * rho_of_L(c * Q / (m * Z), c)
|
|
|
|
|
|
def court2(ch2, QA, QB, JA, JB, t, tau_ref):
|
|
"""To klasser. prio=True: A har forrang (full kapasitet), B får resten.
|
|
prio=False (FIFO): felles kø, lik ventetid. Returnerer R_A, R_B, tau_A, tau_B."""
|
|
Z = ch2["Zfun"](t) if ch2["Zfun"] else ch2["Z"]
|
|
m, c = ch2["m"], ch2["c"]
|
|
if ch2["prio"]:
|
|
RA = capfun(Z, QA, m, c)
|
|
ZB = max(Z - RA, 0.0)
|
|
RB = capfun(ZB, QB, m, c)
|
|
else:
|
|
Qt = QA + QB
|
|
R = capfun(Z, Qt, m, c)
|
|
RA, RB = (R * QA / Qt, R * QB / Qt) if Qt > 0 else (0.0, 0.0)
|
|
tauA = QA / RA if RA > 0 else m
|
|
tauB = QB / RB if RB > 0 else m
|
|
return RA, RB, tauA, tauB
|
|
|
|
|
|
def gfun(tau, tau_ref, tauh):
|
|
D = max(0.0, tau - tau_ref)
|
|
return D / (D + tauh) if D > 0 else 0.0
|
|
|
|
|
|
CH2 = dict(Z=ZD0, m=MD, c=1, prio=True, piA=0.4, thetaA=0.8, thetaB=0.3,
|
|
chi=0.15, mu4=0.12, tauh=9.0, varsigma=VARSIGMA,
|
|
xiA=1.0, xiB=1.0, Zfun=None, xiBfun=None)
|
|
|
|
|
|
def calibrate_m(ch2, p, tau_ref=TAU_REF):
|
|
"""Finn m slik at stasjonær saksbehandlingstid i utgangslikevekten er tau_ref
|
|
(én felles kø, utnyttelse ups_0). For c=1 gir dette MD fra Del 2."""
|
|
Z = ch2["Z"]; J0 = p["gamma"] * p["Lam0"] * p["T0"]
|
|
|
|
def tau_stat(m):
|
|
f = lambda Q: capfun(Z, Q, m, ch2["c"]) + ch2["varsigma"] * Q - J0
|
|
Q0 = brentq(f, 1e-9, 1e4)
|
|
return Q0 / capfun(Z, Q0, m, ch2["c"]) - tau_ref
|
|
return brentq(tau_stat, 1e-4, tau_ref - 1e-6)
|
|
|
|
|
|
def simulate_chain2(tau, ch2, p=P, kp=3.0):
|
|
"""Regime A-styring (P-regulator) med tokjassekø. K_A og K_B er alvorlighets-
|
|
vektet belastning i hver klasse; K = K_A + K_B."""
|
|
n = int(TEND / DT) + 1
|
|
t = np.arange(n) * DT
|
|
piA = ch2["piA"]
|
|
KA = np.full(n, piA); KB = np.full(n, 1.0 - piA)
|
|
U = np.full(n, p["U0"]); L = np.full(n, p["Lam0"]); T = np.full(n, p["T0"])
|
|
QA = np.zeros(n); QB = np.zeros(n); tA = np.zeros(n); tB = np.zeros(n)
|
|
RAs = np.zeros(n); RBs = np.zeros(n)
|
|
# stasjonære køer i utgangslikevekten (felles kø => lik tid; prioritet: A kort)
|
|
J0 = p["gamma"] * p["Lam0"] * p["T0"]
|
|
JA0, JB0 = ch2["xiA"] * piA * J0, ch2["xiB"] * (1 - piA) * J0
|
|
|
|
if ch2["prio"]:
|
|
Qa = brentq(lambda Q: capfun(ch2["Z"], Q, ch2["m"], ch2["c"]) + ch2["varsigma"] * Q - JA0, 1e-12, 1e4)
|
|
ZB0 = ch2["Z"] - capfun(ch2["Z"], Qa, ch2["m"], ch2["c"])
|
|
Qb = brentq(lambda Q: capfun(ZB0, Q, ch2["m"], ch2["c"]) + ch2["varsigma"] * Q - JB0, 1e-12, 1e4)
|
|
else:
|
|
Qt = brentq(lambda Q: capfun(ch2["Z"], Q, ch2["m"], ch2["c"]) + ch2["varsigma"] * Q - (JA0 + JB0), 1e-12, 1e4)
|
|
Qa, Qb = Qt * JA0 / (JA0 + JB0), Qt * JB0 / (JA0 + JB0)
|
|
QA[0], QB[0] = Qa, Qb
|
|
RA0, RB0, trA, trB = court2(ch2, Qa, Qb, JA0, JB0, 0.0, TAU_REF)
|
|
tau_refA, tau_refB = trA, trB # klassevise referansetider (utgangslikevekten)
|
|
lag = int(round(tau / DT))
|
|
for k in range(n - 1):
|
|
K = KA[k] + KB[k]
|
|
Kobs = (KA[k - lag] + KB[k - lag]) if k - lag >= 0 else 1.0
|
|
T[k] = np.clip(p["T0"] + kp * (Kobs - p["Kstar"]), 0.0, p["Tmax"])
|
|
xiB = ch2["xiBfun"](t[k]) if ch2["xiBfun"] else ch2["xiB"]
|
|
h = removal(K, L[k], T[k], p) / K # fjerning per enhet
|
|
if ch2.get("dbr"):
|
|
# trommel-buffer-tau via T: J_A + J_B <= dbr * Z (J er lineær i T)
|
|
Zt = ch2["Zfun"](t[k]) if ch2["Zfun"] else ch2["Z"]
|
|
Jt = h * (ch2["xiA"] * KA[k] + xiB * KB[k])
|
|
if Jt > ch2["dbr"] * Zt:
|
|
T[k] *= ch2["dbr"] * Zt / Jt
|
|
h = removal(K, L[k], T[k], p) / K
|
|
FA, FB = h * KA[k], h * KB[k]
|
|
JA, JB = ch2["xiA"] * FA, xiB * FB
|
|
RA, RB, tA[k], tB[k] = court2(ch2, QA[k], QB[k], JA, JB, t[k], TAU_REF)
|
|
RAs[k], RBs[k] = RA, RB
|
|
gA, gB = gfun(tA[k], tau_refA, ch2["tauh"]), gfun(tB[k], tau_refB, ch2["tauh"])
|
|
gbar = piA * gA + (1 - piA) * gB
|
|
Om = Omega(t[k], p)
|
|
dKA = piA * Om * (1 + ch2["chi"] * gA) - (FA - ch2["thetaA"] * JA) - ch2["thetaA"] * RA - p["delta"] * KA[k]
|
|
dKB = (1 - piA) * Om * (1 + ch2["chi"] * gB) - (FB - ch2["thetaB"] * JB) - ch2["thetaB"] * RB - p["delta"] * KB[k]
|
|
Ueq = np.clip(p["U0"] - p["psiU"] * (K - 1.0), 0.0, 1.0)
|
|
dU = p["kU"] * (Ueq - U[k])
|
|
over = max(0.0, T[k] - p["Tprop"] * K)
|
|
dL = p["kL"] * (p["Lam0"] + p["mu1"] * (U[k] - p["U0"]) - ch2["mu4"] * gbar - L[k]) - p["mu2"] * L[k] * over ** 2
|
|
KA[k + 1] = KA[k] + DT * dKA; KB[k + 1] = KB[k] + DT * dKB
|
|
U[k + 1] = U[k] + DT * dU
|
|
L[k + 1] = np.clip(L[k] + DT * dL, 1e-4, 1 - 1e-4)
|
|
QA[k + 1] = max(QA[k] + DT * (JA - RA - ch2["varsigma"] * QA[k]), 0.0)
|
|
QB[k + 1] = max(QB[k] + DT * (JB - RB - ch2["varsigma"] * QB[k]), 0.0)
|
|
for arr in (T, tA, tB, RAs, RBs):
|
|
arr[-1] = arr[-2]
|
|
return dict(t=t, K=KA + KB, KA=KA, KB=KB, U=U, L=L, T=T, QA=QA, QB=QB, tauA=tA, tauB=tB,
|
|
RA=RAs, RB=RBs, tau_refA=tau_refA, tau_refB=tau_refB)
|
|
|
|
|
|
def summary2(name, r):
|
|
print(f" {name:36s} Kend={r['K'][-1]:.3f} KA_end={r['KA'][-1]:.3f} KB_end={r['KB'][-1]:.3f} "
|
|
f"Lend={r['L'][-1]:.3f} Tend={r['T'][-1]:.2f} tauA(end)={r['tauA'][-1]:5.1f} "
|
|
f"tauB(end)={r['tauB'][-1]:5.1f} tauAmax={r['tauA'].max():4.1f} tauBmax={r['tauB'].max():5.1f} "
|
|
f"QA={r['QA'][-1]:.2f} QB={r['QB'][-1]:.2f}")
|
|
|
|
|
|
# --- (c) rekruttering / smitte ------------------------------------------------
|
|
SOC = dict(eps=0.30, Kh=2.0, pexp=4)
|
|
|
|
|
|
def hill(K, s=SOC):
|
|
return K ** s["pexp"] / (s["Kh"] ** s["pexp"] + K ** s["pexp"])
|
|
|
|
|
|
def hill_d(K, s=SOC):
|
|
pe, Kh = s["pexp"], s["Kh"]
|
|
return pe * K ** (pe - 1) * Kh ** pe / (Kh ** pe + K ** pe) ** 2
|
|
|
|
|
|
def simulate_soc(controller, tau, p, s=SOC, kp=3.0, shock=1.6, t_on=12.0, t_off=36.0):
|
|
"""Aggregert modell med rekrutteringsledd og et FORBIGÅENDE sjokk i Omega
|
|
(faktor shock i [t_on, t_off)). Omega0 rekalibreres slik at K=1 forblir
|
|
likevekt ved T=1. controller: 'P' eller 'fast' (T = T0)."""
|
|
Om0 = p["Omega0"] - s["eps"] * hill(1.0, s)
|
|
n = int(TEND / DT) + 1
|
|
t = np.arange(n) * DT
|
|
K = np.ones(n); U = np.full(n, p["U0"]); L = np.full(n, p["Lam0"]); T = np.full(n, p["T0"])
|
|
lag = int(round(tau / DT))
|
|
for k in range(n - 1):
|
|
Kobs = K[k - lag] if k - lag >= 0 else 1.0
|
|
if controller == "P":
|
|
T[k] = np.clip(p["T0"] + kp * (Kobs - p["Kstar"]), 0.0, p["Tmax"])
|
|
else:
|
|
T[k] = p["T0"]
|
|
Om = Om0 * (shock if t_on <= t[k] < t_off else 1.0) + s["eps"] * hill(K[k], s)
|
|
dK = Om - removal(K[k], L[k], T[k], p) - p["delta"] * K[k]
|
|
Ueq = np.clip(p["U0"] - p["psiU"] * (K[k] - 1.0), 0.0, 1.0)
|
|
dU = p["kU"] * (Ueq - U[k])
|
|
over = max(0.0, T[k] - p["Tprop"] * K[k])
|
|
dL = p["kL"] * (p["Lam0"] + p["mu1"] * (U[k] - p["U0"]) - L[k]) - p["mu2"] * L[k] * over ** 2
|
|
K[k + 1] = K[k] + DT * dK
|
|
U[k + 1] = U[k] + DT * dU
|
|
L[k + 1] = np.clip(L[k] + DT * dL, 1e-4, 1 - 1e-4)
|
|
T[-1] = T[-2]
|
|
return t, K, U, L, T
|
|
|
|
|
|
def soc_equilibria(p, s=SOC, T=1.0, L=0.8, Tmax=None, kp=None):
|
|
"""Likevekter for dK/dt = 0 ved fast T (eller P-regel hvis kp gis), Lambda = L."""
|
|
Om0 = p["Omega0"] - s["eps"] * hill(1.0, s)
|
|
|
|
def Tof(K):
|
|
if kp is None:
|
|
return T
|
|
return np.clip(p["T0"] + kp * (K - p["Kstar"]), 0.0, Tmax if Tmax else p["Tmax"])
|
|
|
|
def f(K):
|
|
return Om0 + s["eps"] * hill(K, s) - removal(K, L, Tof(K), p) - p["delta"] * K
|
|
Ks = np.linspace(0.01, 12.0, 24000); fv = np.array([f(x) for x in Ks])
|
|
return [brentq(f, Ks[i], Ks[i + 1]) for i in range(len(Ks) - 1) if fv[i] * fv[i + 1] < 0]
|
|
|
|
|
|
# =============================================================================
|
|
# Del 4: endogene tiltak -- fra observasjon til beslutning (avsnitt 8.5)
|
|
# =============================================================================
|
|
# Tiltaksvektoren v = (T, xi_D, Z_D, q) bestemmes av beslutningsregelen Psi(x_hat, O), lign. (9).
|
|
# Observasjoner O: K med forsinkelse tau (politiets sløyfe, som før); saksbehandlingstid,
|
|
# innstrømning, kø og kapasitet i domstolen med forsinkelse lagD; tillit fra årlig
|
|
# undersøkelse publisert lag_survey måneder etter årsskiftet.
|
|
# - "sched": faste tidsplaner for xi og Z (gjenskaper scenario 2-5 eksakt)
|
|
# - "desentral": påtalemyndigheten justerer xi etter observert saksbehandlingstid
|
|
# (integralregel med dødsone); budsjettet bestiller kapasitet etter
|
|
# observert innstrømning siste 12 mnd (bestilt kapasitet regnes med);
|
|
# ingen av dem kjenner den andres planer
|
|
# - "desentral_prog": som over, men budsjettet dimensjonerer etter prognosen for
|
|
# strømmen i gjenopprettet likevekt (uten å kjenne xi-planen)
|
|
# - "koord": felles kapasitetsbalanse når domstolen er observert overbelastet
|
|
# (ups_D > 1): xi settes slik at innstrømning pluss
|
|
# restanseavvikling holdes innenfor ups* Z (trommel-buffer-tau på
|
|
# påtalebeslutningene); budsjettet bestiller bare det
|
|
# underordningen ikke tar unna, dimensjonert etter prognosen
|
|
# Kapasitet realiseres lineært i [t+L1, t+L2] etter vedtak; nedbygging skjer ved naturlig
|
|
# avgang med tidskonstant tau_down og bare når kapasiteten er overflødig (dødsone).
|
|
ENDO = dict(
|
|
lagD=3.0, # observasjonsforsinkelse for kjedestatistikk (mnd)
|
|
xi_min=0.8, # gulv: høyst 20 % av sakene til forenklede reaksjoner
|
|
tau_thr=6.0, # terskel for observert saksbehandlingstid (desentral xi-regel)
|
|
k_xi=0.2 / 12, # integralforsterkning: 0,2 per år når tiden er dobbelt så lang som terskelen
|
|
tau_xi=6.0, # tidskonstant for endring av påtalepraksis
|
|
ups_star=0.95, # målutnyttelse ved dimensjonering
|
|
ups_hi=1.0, ups_lo=0.85, # dødsone: bestill ved observert overbelastning, bygg ned under 0,85
|
|
kZ=1.0, # budsjettforsterkning (1 = hele gapet bestilles)
|
|
L1=6.0, L2=24.0, # realisering av bestilt kapasitet: lineært i [t+L1, t+L2]
|
|
tau_down=36.0, # tidskonstant for nedbygging (naturlig avgang)
|
|
T_clear=24.0, # koordinert regel: restansen skal avvikles over denne tiden
|
|
H_pred=36.0, # prognosehorisont (mnd)
|
|
period=12.0, # budsjettsyklus
|
|
qloop=False, mu3=0.5, Lam_star=0.78, k_q=1.0 / 12, q_max=0.2, lag_survey=3.0,
|
|
)
|
|
|
|
|
|
def deriv_chain_xz(K, U, L, T, Q, Z, xi, q, t, p, ch, tau_ref, Om0=None, mu3=0.0):
|
|
"""Som deriv_chain, men med Z, xi og q som eksplisitte argumenter."""
|
|
F = removal(K, L, T, p)
|
|
JD = xi * F
|
|
RD, tauD, dQ, g = court(dict(ch, Z=Z, Zfun=None), Q, JD, t, tau_ref)
|
|
Om = (Omega(t, p) if Om0 is None else Om0) * (1.0 + ch["chi"] * g)
|
|
dK = Om - (1.0 - xi * ch["theta"]) * F - ch["theta"] * RD - p["delta"] * K
|
|
Ueq = np.clip(p["U0"] - p["psiU"] * (K - 1.0), 0.0, 1.0)
|
|
dU = p["kU"] * (Ueq - U)
|
|
over = max(0.0, T - p["Tprop"] * K)
|
|
dL = (p["kL"] * (p["Lam0"] + p["mu1"] * (U - p["U0"]) + mu3 * q - ch["mu4"] * g - L)
|
|
- p["mu2"] * L * over ** 2)
|
|
return dK, dU, dL, dQ, tauD, JD, RD, F
|
|
|
|
|
|
def predict_F(K0, U0, L0, Q0, Om0, p, ch, tau_ref, H, kp, lag_steps):
|
|
"""Prognose: kjør den interne modellen uten flaskehals (Z -> inf, alle saker til
|
|
domstol) H måneder fram fra observert tilstand, med samme styringsregel for T;
|
|
returner fjerningsstrømmen F i den gjenopprettede likevekten."""
|
|
K, U, L, Q = K0, U0, L0, Q0
|
|
Kq = [K0] * (lag_steps + 1)
|
|
F = 0.0
|
|
for j in range(int(H / DT)):
|
|
Kobs = Kq[-lag_steps - 1]
|
|
T = np.clip(p["T0"] + kp * (Kobs - p["Kstar"]), 0.0, p["Tmax"])
|
|
dK, dU, dL, dQ, _, _, _, F = deriv_chain_xz(K, U, L, T, Q, np.inf, 1.0, 0.0, 0.0, p, ch, tau_ref, Om0=Om0)
|
|
K += DT * dK; U += DT * dU; L = np.clip(L + DT * dL, 1e-4, 1 - 1e-4); Q = max(Q + DT * dQ, 0.0)
|
|
Kq.append(K)
|
|
return F
|
|
|
|
|
|
def simulate_endo(policy, tau=3.0, p=P, ch=CH0, e=None, kp=3.0, tend=TEND, xi_sched=None, Z_sched=None):
|
|
e = dict(ENDO, **(e or {}))
|
|
n = int(tend / DT) + 1
|
|
t = np.arange(n) * DT
|
|
K = np.ones(n); U = np.full(n, p["U0"]); L = np.full(n, p["Lam0"]); T = np.full(n, p["T0"])
|
|
Q = np.zeros(n); tauD = np.zeros(n); JD = np.zeros(n); RD = np.zeros(n); Fa = np.zeros(n)
|
|
Z = np.full(n, ch["Z"]); xi = np.ones(n); q = np.zeros(n)
|
|
Q0, tau_ref = init_chain(dict(ch, Zfun=None, xifun=None), p)
|
|
Q[0] = Q0
|
|
lag = int(round(tau / DT)); lagD = int(round(e["lagD"] / DT))
|
|
mu3 = e["mu3"] if e["qloop"] else 0.0
|
|
xi_tgt, xi_int, q_tgt = 1.0, 1.0, 0.0
|
|
binding = False # koord: er domstolen identifisert som bindende skranke (ToC trinn 1)?
|
|
orders = [] # bestilt kapasitet: (t_start, t_slutt, dZ), realiseres lineært
|
|
Z_down = None # nedbyggingsmål
|
|
log = [] # budsjettvedtak
|
|
for k in range(n - 1):
|
|
tk = t[k]
|
|
Kobs = K[k - lag] if k - lag >= 0 else 1.0
|
|
T[k] = np.clip(p["T0"] + kp * (Kobs - p["Kstar"]), 0.0, p["Tmax"])
|
|
i = k - lagD # forsinkede observasjoner av kjeden
|
|
if i >= 0:
|
|
tau_hat, J_hat, Q_hat, F_hat = tauD[i], JD[i], Q[i], Fa[i]
|
|
else:
|
|
tau_hat, J_hat, Q_hat, F_hat = tau_ref, 0.1, Q0, 0.1
|
|
at_month = abs(tk - round(tk)) < 1e-9
|
|
if policy == "sched":
|
|
xi[k] = xi_sched(tk) if xi_sched else 1.0
|
|
Z[k] = Z_sched(tk) if Z_sched else ch["Z"]
|
|
else:
|
|
Q_ref = (tau_ref - ch["m"]) * Z[k] # kø som svarer til referansetiden
|
|
Q_exc = max(Q_hat - Q_ref, 0.0)
|
|
# ---- påtalemyndighetens sløyfe (månedlig) ----
|
|
if at_month and k > 0:
|
|
if policy in ("desentral", "desentral_prog"):
|
|
if tau_hat > e["tau_thr"]:
|
|
xi_int -= e["k_xi"] * min(tau_hat / e["tau_thr"] - 1.0, 3.0)
|
|
elif tau_hat < tau_ref:
|
|
xi_int += e["k_xi"] * min(1.0 - tau_hat / tau_ref, 3.0)
|
|
xi_int = float(np.clip(xi_int, e["xi_min"], 1.0))
|
|
xi_tgt = xi_int
|
|
else: # koord: trommel-buffer-tau på påtalebeslutningene når skranken binder
|
|
binding = binding or F_hat > e["ups_hi"] * Z[k]
|
|
xi_tgt = float(np.clip((e["ups_star"] * Z[k] - Q_exc / e["T_clear"]) / max(F_hat, 1e-9),
|
|
e["xi_min"], 1.0)) if binding else 1.0
|
|
if e["qloop"]: # kvalitetssløyfen: årlig undersøkelse, publisert med forsinkelse
|
|
ks = int(round(np.floor((tk - e["lag_survey"]) / 12.0) * 12.0 / DT))
|
|
L_hat = L[ks] if ks >= 0 else p["Lam0"]
|
|
q_tgt = float(np.clip(q_tgt + e["k_q"] * (e["Lam_star"] - L_hat), 0.0, e["q_max"]))
|
|
# ---- budsjettsløyfen (årlig) ----
|
|
if at_month and k > 0 and abs(tk % e["period"]) < 1e-9:
|
|
pending = sum(dZ * (1 - np.clip((tk - a) / (b - a), 0, 1)) for a, b, dZ in orders)
|
|
Z_comm = Z[k] + pending # realisert pluss bevilget kapasitet
|
|
if policy == "desentral":
|
|
i0 = max(i - int(12 / DT), 0)
|
|
Jbar = JD[i0:i + 1].mean() if i > 0 else J_hat
|
|
need, surplus_ref = Jbar, Jbar
|
|
sig = f"ups_obs={Jbar / Z[k]:.2f}"
|
|
else:
|
|
i_ = max(i, 0)
|
|
F1 = predict_F(K[i_], U[i_], L[i_], Q[i_], Omega(t[i_], p), p, ch, tau_ref, e["H_pred"], kp, lag)
|
|
xi_pl = e["xi_min"] if policy == "koord" else 1.0
|
|
need, surplus_ref = xi_pl * F1 + Q_exc / e["T_clear"], F1
|
|
sig = f"F_pred={F1:.3f} Q_exc={Q_exc:.2f}"
|
|
if need > e["ups_hi"] * Z_comm:
|
|
dZ = e["kZ"] * (need / e["ups_star"] - Z_comm)
|
|
orders.append((tk + e["L1"], tk + e["L2"], dZ)); Z_down = None
|
|
log.append(f"mnd {tk:.0f}: bestill +{dZ / ch['Z'] * 100:.0f} % ({sig})")
|
|
elif surplus_ref < e["ups_lo"] * Z[k] and pending < 1e-9:
|
|
Z_down = surplus_ref / e["ups_star"]
|
|
log.append(f"mnd {tk:.0f}: bygg ned mot {Z_down / ch['Z']:.2f} ({sig})")
|
|
# ---- realisering ----
|
|
rate = sum(dZ / (b - a) for a, b, dZ in orders if a <= tk < b)
|
|
if Z_down is not None and rate == 0.0 and Z[k] > Z_down:
|
|
rate -= (Z[k] - Z_down) / e["tau_down"]
|
|
Z[k + 1] = Z[k] + DT * rate
|
|
xi[k + 1] = xi[k] + DT * (xi_tgt - xi[k]) / e["tau_xi"]
|
|
q[k + 1] = q[k] + DT * (q_tgt - q[k]) / 12.0
|
|
dK, dU, dL, dQ, tauD[k], JD[k], RD[k], Fa[k] = deriv_chain_xz(
|
|
K[k], U[k], L[k], T[k], Q[k], Z[k], xi[k], q[k], tk, p, ch, tau_ref, mu3=mu3)
|
|
K[k + 1] = K[k] + DT * dK
|
|
U[k + 1] = U[k] + DT * dU
|
|
L[k + 1] = np.clip(L[k] + DT * dL, 1e-4, 1 - 1e-4)
|
|
Q[k + 1] = max(Q[k] + DT * dQ, 0.0)
|
|
if policy == "sched":
|
|
xi[-1] = xi_sched(t[-1]) if xi_sched else 1.0
|
|
Z[-1] = Z_sched(t[-1]) if Z_sched else ch["Z"]
|
|
for arr in (T, tauD, JD, RD, Fa):
|
|
arr[-1] = arr[-2]
|
|
return dict(t=t, K=K, U=U, L=L, T=T, Q=Q, tauD=tauD, JD=JD, RD=RD, Z=Z, xi=xi, q=q, tau_ref=tau_ref, log=log)
|
|
|
|
|
|
def summary_endo(name, r, Z0=ZD0):
|
|
t = r["t"]; m = t >= 12
|
|
frafall = np.trapezoid(VARSIGMA * r["Q"], t) # kumulativt frafall (normaliserte saker)
|
|
capyears = np.trapezoid(np.maximum(r["Z"] - Z0, 0) / Z0, t) / 12 # kapasitetsår over referansen
|
|
print(f" {name:36s} Kmean={r['K'][m].mean():.3f} Kend={r['K'][-1]:.3f} Lmin={r['L'].min():.3f} "
|
|
f"Lend={r['L'][-1]:.3f} tauDmax={r['tauD'].max():5.1f} tauDend={r['tauD'][-1]:4.1f} "
|
|
f"Zmax={r['Z'].max() / Z0:.2f} Zend={r['Z'][-1] / Z0:.2f} xi_min={r['xi'].min():.2f} "
|
|
f"xi_end={r['xi'][-1]:.2f} frafall={frafall:.2f} kap.år={capyears:.2f} Tend={r['T'][-1]:.2f}")
|
|
return dict(Kmean=r['K'][m].mean(), frafall=frafall, capyears=capyears)
|
|
|
|
|
|
def pure_capacity_loop(kZ, L1=6.0, L2=24.0, supply_line=True, tend=600.0, J=0.134, ups=0.95):
|
|
"""Ren kapasitetssløyfe med gitt innstrømning J: årlig vedtak kZ*(J/ups - Z_comm),
|
|
realisert lineært i [t+L1, t+L2]. supply_line=False: bevilget, ikke realisert
|
|
kapasitet regnes ikke med (Z_comm = Z). Returnerer Z/Z_behov over tid."""
|
|
n = int(tend / DT) + 1
|
|
t = np.arange(n) * DT; Z = np.full(n, ZD0); orders = []
|
|
for k in range(n - 1):
|
|
tk = t[k]
|
|
if abs(tk - round(tk)) < 1e-9 and tk > 0 and tk % 12 == 0:
|
|
pending = sum(dZ * (1 - np.clip((tk - a) / (b - a), 0, 1)) for a, b, dZ in orders) if supply_line else 0.0
|
|
orders.append((tk + L1, tk + L2, kZ * (J / ups - (Z[k] + pending))))
|
|
Z[k + 1] = Z[k] + DT * sum(dZ / (b - a) for a, b, dZ in orders if a <= tk < b)
|
|
return t, Z / (J / ups)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
import os
|
|
os.makedirs("data", exist_ok=True)
|
|
print("=== Del 1: styringsregimer uten kjede (tabell 3, avsnitt 8.1 og 9.3) ===")
|
|
a = P["gamma"] * P["Lam0"] * P["T0"] + P["delta"]
|
|
for kp in (1.0, 3.0, 6.0):
|
|
b = P["gamma"] * P["Lam0"] * 1.0 * kp
|
|
tc = critical_delay(a, b)
|
|
per = 2 * np.pi / np.sqrt(b ** 2 - a ** 2) if b > a else np.inf
|
|
print(f"kp={kp}: a={a:.3f} b={b:.3f} tau_c={tc:.2f} mnd periode={per:.1f} mnd")
|
|
c = P["gamma"] * P["Lam0"]
|
|
print(f"stabil for alle tau hvis kp <= a/c = {a / c:.2f}")
|
|
tau = 11.0
|
|
kp_max = brentq(lambda k: critical_delay(a, c * k) - tau, a / c + 1e-6, 50.0)
|
|
print(f"tau={tau}: største stabile kp = {kp_max:.2f}; tommelfingerregel kp < {np.pi / (2 * c * tau):.2f}")
|
|
|
|
runs = {"p3": ("P", 3.0), "p6": ("P", 6.0), "p11": ("P", 11.0), "mpc11": ("MPC", 11.0)}
|
|
base = {}
|
|
for name, (ctrl, tau) in runs.items():
|
|
t, K, U, L, T = simulate(ctrl, tau)
|
|
base[name] = (t, K, U, L, T)
|
|
write_csv(f"data/{name}.dat", [t, K, U, L, T], "t K U L T")
|
|
late = t > 60
|
|
print(f"{name:6s} tau={tau:4.1f}: Kmax={K.max():.3f} Kmin={K[t>12].min():.3f} "
|
|
f"Tmin={T.min():.2f} Tmax={T.max():.2f} Lmin={L.min():.3f} "
|
|
f"amp(K,t>60)={K[late].max()-K[late].min():.3f} Kend={K[-1]:.3f} Lend={L[-1]:.3f}")
|
|
|
|
# regime E: avdempet proporsjonal regel (kp = 1,2) med samme forsinkelse som C
|
|
t, K, U, L, T = simulate("P", 11.0, kp=1.2)
|
|
write_csv("data/p11d.dat", [t, K, U, L, T], "t K U L T")
|
|
late = t > 60
|
|
print(f"p11d tau=11.0 kp=1.2: Kmax={K.max():.3f} Tmax={T.max():.2f} Lmin={L.min():.3f} "
|
|
f"amp(K,t>60)={K[late].max()-K[late].min():.3f} Kend={K[-1]:.3f} Lend={L[-1]:.3f}")
|
|
|
|
# stabilitetskart: kritisk forsinkelse som funksjon av forsterkning kp
|
|
kps = np.linspace(0.05, 12, 400)
|
|
rows = []
|
|
for kp in kps:
|
|
b = P["gamma"] * P["Lam0"] * kp
|
|
tc = critical_delay(a, b)
|
|
rows.append((kp, min(tc, 60.0)))
|
|
np.savetxt("data/stab.dat", np.array(rows), header="kp tauc", comments="", fmt="%.4f")
|
|
# tommelfingerregel b*tau < pi/2 (a -> 0)
|
|
rows2 = [(kp, min(np.pi / 2 / (P["gamma"] * P["Lam0"] * kp), 60.0)) for kp in kps]
|
|
np.savetxt("data/thumb.dat", np.array(rows2), header="kp tauc", comments="", fmt="%.4f")
|
|
|
|
# ------------------------------------------------------------------------
|
|
print("\n=== Del 2: straffesakskjeden koblet til dynamikken (avsnitt 7.3, 7.6, 8.2) ===")
|
|
r0 = simulate_chain("P", 3.0, CH0)
|
|
Q0 = r0["Q"][0]
|
|
print(f"Z_D={ZD0:.4f} (ups_0={UPS0}), m_D={MD:.3f} mnd (kalibrert), tau_D^ref={r0['tau_ref']:.2f} mnd, "
|
|
f"Q_D(0)={Q0:.3f}, frafall før sjokk={VARSIGMA*Q0/J0*100:.1f} % av J_D, "
|
|
f"faktisk utnyttelse R_D/Z_D={(J0-VARSIGMA*Q0)/ZD0:.3f}")
|
|
ups_c = (J0 - VARSIGMA * Q0) / ZD0
|
|
print(f"lineær relaksasjonstid for køen: 1/((1-ups)^2/m_D + varsigma) = "
|
|
f"{1/((1-ups_c)**2/MD + VARSIGMA):.0f} mnd (uten frafall: {MD/(1-ups_c)**2:.0f} mnd)")
|
|
print(f"M/M/1-sammenheng: tau_D = m_D/(1-ups): ups=0.8 -> {MD/0.2:.2f} mnd, ups=0.95 -> {MD/0.05:.2f} mnd")
|
|
|
|
# --- grensetilfelle-kontroll ---------------------------------------------
|
|
print("\n--- grensetilfelle: reproduksjon av regime A-D ---")
|
|
for name, (ctrl, tau) in runs.items():
|
|
t, K, U, L, T = base[name]
|
|
r = simulate_chain(ctrl, tau, dict(CH0, Z=np.inf), mpc_knows_chain=False)
|
|
print(f" Z_D=inf, {name:6s}: max|dK|={np.abs(r['K']-K).max():.2e} max|dL|={np.abs(r['L']-L).max():.2e} "
|
|
f"max|dT|={np.abs(r['T']-T).max():.2e}")
|
|
t, K, U, L, T = base["p3"]
|
|
r = simulate_chain("P", 3.0, dict(CH0, theta=0.0, chi=0.0, mu4=0.0))
|
|
print(f" Z_D=Z_D0, tilbakekoblinger av (theta=chi=mu4=0): max|dK|={np.abs(r['K']-K).max():.2e} "
|
|
f"max|dL|={np.abs(r['L']-L).max():.2e} (køen beregnes, men virker ikke tilbake)")
|
|
for Zbig in (1.0, 10.0):
|
|
r = simulate_chain("P", 3.0, dict(CH0, Z=Zbig))
|
|
print(f" Z_D={Zbig:4.1f} (endelig, m_D={MD:.3f}): max|dK|={np.abs(r['K']-K).max():.1e} "
|
|
f"max|dL|={np.abs(r['L']-L).max():.1e} tau_D^ref={r['tau_ref']:.3f}")
|
|
|
|
# --- scenarioer (regime A: kp=3, tau=3) ---------------------------------
|
|
T_TILTAK = 24.0 # beslutning 12 måneder etter sjokket (restansen er observert)
|
|
T_UTV = 24.0 # kapasitetsøkningen realiseres lineært over 24 måneder
|
|
XI_UND = 0.8 # underordning: 20 % av sakene til forenklede reaksjoner
|
|
ZD_UTV = 1.34 # utvidelse: +34 % kapasitet, til 0.141 (jf. regneeksempelet i avsnitt 7.6)
|
|
xi_und = lambda t: XI_UND if t >= T_TILTAK else 1.0
|
|
Z_utv = lambda t: ZD0 * (1 + (ZD_UTV - 1) * np.clip((t - T_TILTAK) / T_UTV, 0, 1))
|
|
scen = {
|
|
"s1": ("1 ingen flaskehals (= regime A)", dict(CH0, Z=np.inf)),
|
|
"s2": ("2 flaskehals, ingen tiltak", CH0),
|
|
"s3": ("3 underordne (xi_D 1->0.8 fra mnd 24)", dict(CH0, xifun=xi_und)),
|
|
"s4": ("4 utvide (Z_D +34 % over mnd 24-48)", dict(CH0, Zfun=Z_utv)),
|
|
"s5": ("5 underordne + utvide", dict(CH0, Zfun=Z_utv, xifun=xi_und)),
|
|
}
|
|
print("\n--- scenarioer, regime A (kp=3, tau=3), fig:kjedesim / tabell ---")
|
|
res = {}
|
|
for key, (label, ch) in scen.items():
|
|
r = simulate_chain("P", 3.0, ch)
|
|
res[key] = r
|
|
write_csv(f"data/kjede_{key}.dat", [r["t"], r["K"], r["U"], r["L"], r["T"], r["Q"], r["tauD"]],
|
|
"t K U L T Q tauD")
|
|
summary(label, r)
|
|
i24, i48 = int(24 / DT), int(48 / DT)
|
|
for key, r in res.items():
|
|
print(f" {key}: tauD@24={r['tauD'][i24]:.1f} tauD@48={r['tauD'][i48]:.1f} K@48={r['K'][i48]:.3f} "
|
|
f"JD_end={r['JD'][-1]:.3f} RD_end={r['RD'][-1]:.3f} frafall_end={VARSIGMA*r['Q'][-1]:.3f} "
|
|
f"Umin={r['U'].min():.3f}")
|
|
r1 = res["s1"]
|
|
print(f" s1: J_D i ny likevekt={r1['JD'][-1]:.3f} (ups_D={r1['JD'][-1]/ZD0:.2f}); med xi=0.8: "
|
|
f"{0.8*r1['JD'][-1]:.3f} (ups_D={0.8*r1['JD'][-1]/ZD0:.2f}); Z_D for ups=0.95: {r1['JD'][-1]/0.95:.3f} "
|
|
f"(+{(r1['JD'][-1]/0.95/ZD0-1)*100:.0f} %), med xi=0.8: {0.8*r1['JD'][-1]/0.95:.3f} (+{(0.8*r1['JD'][-1]/0.95/ZD0-1)*100:.0f} %)")
|
|
|
|
# --- trommel-buffer-tau og MPC --------------------------------------------
|
|
print("\n--- trommel-buffer-tau (T begrenses slik at J_D <= 0.95 Z_D) og MPC (tau=3) ---")
|
|
summary("6 TBT, xi_D=1", simulate_chain("P", 3.0, dict(CH0, dbr=0.95)))
|
|
summary("7 TBT + underordne", simulate_chain("P", 3.0, dict(CH0, dbr=0.95, xifun=xi_und)))
|
|
summary("MPC uten kjedemodell (flaskehals)", simulate_chain("MPC", 3.0, CH0, mpc_knows_chain=False))
|
|
summary("MPC med kjedemodell (flaskehals)", simulate_chain("MPC", 3.0, CH0, mpc_knows_chain=True))
|
|
summary("MPC med kjedemodell + underordne", simulate_chain("MPC", 3.0, dict(CH0, xifun=xi_und), mpc_knows_chain=True))
|
|
|
|
# --- følsomhet ------------------------------------------------------------
|
|
print("\n--- følsomhet for celerity-parameteren chi_D (tabell): K_end / Lambda_end / tau_D(end) ---")
|
|
for key in ("s2", "s3", "s5"):
|
|
label, ch = scen[key]
|
|
row = []
|
|
for chi in (0.0, 0.15, 0.30):
|
|
r = simulate_chain("P", 3.0, dict(ch, chi=chi))
|
|
row.append(f"chi={chi:.2f}: {r['K'][-1]:.3f} / {r['L'][-1]:.3f} / {r['tauD'][-1]:.1f}")
|
|
print(f" {label:38s} " + " | ".join(row))
|
|
print("--- følsomhet for vartheta og mu_4 (scenario 2) ---")
|
|
for th in (0.0, 0.25, 0.5, 0.75):
|
|
r = simulate_chain("P", 3.0, dict(CH0, theta=th))
|
|
print(f" vartheta={th:.2f}: Kend={r['K'][-1]:.3f} Lend={r['L'][-1]:.3f} tauD={r['tauD'][-1]:.1f} Tend={r['T'][-1]:.2f}")
|
|
for m4 in (0.0, 0.06, 0.12, 0.24):
|
|
r = simulate_chain("P", 3.0, dict(CH0, mu4=m4))
|
|
print(f" mu4={m4:.2f}: Kend={r['K'][-1]:.3f} Lend={r['L'][-1]:.3f} tauD={r['tauD'][-1]:.1f}")
|
|
for vs in (0.0, 0.01, 0.02):
|
|
r = simulate_chain("P", 3.0, dict(CH0, varsigma=vs))
|
|
print(f" varsigma={vs:.2f}: Kend={r['K'][-1]:.3f} Lend={r['L'][-1]:.3f} tauD={r['tauD'][-1]:.1f}")
|
|
|
|
# --- stabilitet med kø (regime B og C) -------------------------------------
|
|
print("\n--- stabilitet: regime B (tau=6) og C (tau=11) med mettet domstol ---")
|
|
for tau in (6.0, 11.0):
|
|
ro = simulate_chain("P", tau, dict(CH0, Z=np.inf))
|
|
r2 = simulate_chain("P", tau, CH0)
|
|
r5 = simulate_chain("P", tau, dict(CH0, Zfun=Z_utv, xifun=xi_und))
|
|
for lab, r in (("uten kjede", ro), ("flaskehals (scen. 2)", r2), ("underordne+utvide (scen. 5)", r5)):
|
|
mx, mn = extrema(r)
|
|
summary(f"tau={tau:4.1f} {lab}", r)
|
|
print(f" topper {mx} bunner {mn}")
|
|
r2 = res["s2"]
|
|
Lb, Tb, Kb = r2["L"][-1], r2["T"][-1], r2["K"][-1]
|
|
cb = P["gamma"] * Lb
|
|
s = 1 - CH0["xi"] * CH0["theta"]
|
|
a_sat, b_sat = s * cb * Tb + P["delta"], s * cb * Kb * 3.0
|
|
a_un, b_un = cb * Tb + P["delta"], cb * Kb * 3.0
|
|
print(f" lineær analyse i sluttilstanden til scen. 2 (Lambda={Lb:.3f}, T={Tb:.2f}, K={Kb:.2f}): "
|
|
f"mettet domstol a'={a_sat:.3f} b'={b_sat:.3f} tau_c'={critical_delay(a_sat, b_sat):.1f} mnd; "
|
|
f"samme tilstand uten metning a={a_un:.3f} b={b_un:.3f} tau_c={critical_delay(a_un, b_un):.1f} mnd")
|
|
|
|
# ========================================================================
|
|
print("\n=== Del 3: robusthet (avsnitt 8.3 og 8.4) ===")
|
|
# --- (a) politiets kapasitet -------------------------------------------
|
|
print("\n--- (a) politimetning: lineær analyse i referanselikevekten (b uendret, a = 0.10 s + delta, s = 1/(1+1/K_sat)) ---")
|
|
rbar = P["gamma"] * P["Lam0"] * P["T0"] # fjerning per enhet i referansen (0.10)
|
|
for Ks in (1.0, 1.5, 3.0, np.inf):
|
|
s_ = 1.0 / (1.0 + 1.0 / Ks)
|
|
a_s, b_s = rbar * s_ + P["delta"], rbar * 3.0
|
|
kp_all = a_s / rbar
|
|
kp_max = brentq(lambda k: critical_delay(a_s, rbar * k) - 11.0, kp_all + 1e-6, 50.0)
|
|
print(f" K_sat={Ks:4}: s={s_:.3f} a={a_s:.3f} b={b_s:.3f} tau_c={critical_delay(a_s, b_s):.3f} mnd "
|
|
f"(periode {2*np.pi/np.sqrt(b_s**2-a_s**2):.1f}); stabil for alle tau hvis kp<={kp_all:.2f}; "
|
|
f"største stabile kp ved tau=11: {kp_max:.2f}")
|
|
print(" oppklaringsintensitet h(K)/h(1) = (1+1/K_sat)/(1+K/K_sat) ved fast T:")
|
|
for Ks in (1.5, 3.0):
|
|
print(f" K_sat={Ks}: K=1.3 -> {(1+1/Ks)/(1+1.3/Ks):.3f}, K=1.5 -> {(1+1/Ks)/(1+1.5/Ks):.3f}, K=2 -> {(1+1/Ks)/(1+2/Ks):.3f}")
|
|
|
|
print("--- (a) grensetilfelle K_sat=inf: regime A-C og scenario 2 gjenskapes ---")
|
|
for name, (ctrl, tau) in runs.items():
|
|
if ctrl != "P":
|
|
continue
|
|
t, K, U, L, T = base[name]
|
|
t2, K2, U2, L2, T2 = simulate("P", tau, p=dict(P, ksat=np.inf))
|
|
print(f" K_sat=inf {name}: max|dK|={np.abs(K2-K).max():.2e} max|dT|={np.abs(T2-T).max():.2e}")
|
|
r = simulate_chain("P", 3.0, CH0, p=dict(P, ksat=np.inf))
|
|
print(f" K_sat=inf scen. 2: max|dK|={np.abs(r['K']-res['s2']['K']).max():.2e}")
|
|
|
|
print("--- (a) regime A/B/C med politimetning (K_sat = 1.5 og 3; gamma rekalibrert) ---")
|
|
sat_runs = {}
|
|
for Ks in (3.0, 1.5):
|
|
pk = dict(P, ksat=Ks)
|
|
for name, tau in (("A", 3.0), ("B", 6.0), ("C", 11.0)):
|
|
t, K, U, L, T = simulate("P", tau, p=pk)
|
|
sat_runs[(Ks, name)] = (t, K, U, L, T)
|
|
late = t > 60
|
|
Kend, Tend, Lend = K[-1], T[-1], L[-1]
|
|
h_rel = (1 + 1 / Ks) / (1 + K.max() / Ks)
|
|
# lineær analyse i sluttlikevekten (bare meningsfull når amp ~ 0)
|
|
Lb = L[-1]; cb = P["gamma"] * Lb * (1 + 1 / Ks)
|
|
sK = 1 / (1 + Kend / Ks)
|
|
a_e = cb * Tend * sK ** 2 + P["delta"]; b_e = cb * Kend * sK * 3.0
|
|
mx, mn = extrema(dict(K=K, t=t))
|
|
print(f" K_sat={Ks:3} regime {name} tau={tau:4.1f}: Kmax={K.max():.3f} Kmin={K[t>12].min():.3f} "
|
|
f"Tmax={T.max():.2f} Lmin={L.min():.3f} amp(K,t>60)={K[late].max()-K[late].min():.3f} "
|
|
f"Kend={Kend:.3f} Lend={Lend:.3f} h(Kmax)/h(1)={h_rel:.3f} "
|
|
f"| lin. etter sjokk: a={a_e:.3f} b={b_e:.3f} tau_c={critical_delay(a_e,b_e):.1f}")
|
|
if name in ("B", "C"):
|
|
print(f" topper {mx[:6]} bunner {mn[:6]}")
|
|
# figurdata: regime B og C med og uten metning
|
|
t, K, U, L, T = base["p6"]; tb3 = sat_runs[(3.0, "B")]; tb15 = sat_runs[(1.5, "B")]
|
|
write_csv("data/sat_B.dat", [t, K, tb3[1], tb15[1]], "t K Ks3 Ks15")
|
|
t, K, U, L, T = base["p11"]; tc3 = sat_runs[(3.0, "C")]; tc15 = sat_runs[(1.5, "C")]
|
|
write_csv("data/sat_C.dat", [t, K, tc3[1], tc15[1]], "t K Ks3 Ks15")
|
|
|
|
print("--- (a) MPC (regime D, tau=11) når det styrte systemet er mettet, men MPC-ens interne modell ikke er det ---")
|
|
for Ks in (3.0, 1.5):
|
|
t, K, U, L, T = simulate("MPC", 11.0, p=dict(P, ksat=Ks))
|
|
late = t > 60
|
|
print(f" K_sat={Ks:3}: Kmax={K.max():.3f} Kend={K[-1]:.3f} Tmax={T.max():.2f} Tend={T[-1]:.2f} Lmin={L.min():.3f} "
|
|
f"Lend={L[-1]:.3f} amp(K,t>60)={K[late].max()-K[late].min():.3f} (K_sat=inf: Kmax=1.224 Kend=1.031 Lmin=0.783)")
|
|
print("--- (a) kjedescenarioer (regime A) med politimetning: Kend / Lend / tauD(end) / Tend / R_D-andel ---")
|
|
for Ks in (np.inf, 3.0, 1.5):
|
|
pk = dict(P, ksat=Ks)
|
|
row = []
|
|
for key in ("s2", "s3", "s4", "s5"):
|
|
label, ch = scen[key]
|
|
r = simulate_chain("P", 3.0, ch, p=pk)
|
|
row.append(f"{key}: {r['K'][-1]:.3f}/{r['L'][-1]:.3f}/{r['tauD'][-1]:4.1f}/{r['T'][-1]:.2f}")
|
|
print(f" K_sat={Ks:4}: " + " | ".join(row))
|
|
|
|
# --- (b) prioritert tokjassekø og flerserver -----------------------------
|
|
print("\n--- (b) domstol med to sakstyper: kontrolltall ---")
|
|
r_fifo = simulate_chain2(3.0, dict(CH2, prio=False, thetaA=0.5, thetaB=0.5))
|
|
print(f" FIFO, theta_A=theta_B=0.5, c=1: max|dK| mot scenario 2 = {np.abs(r_fifo['K']-res['s2']['K']).max():.2e}, "
|
|
f"max|dtau| = {np.abs(r_fifo['tauA']-res['s2']['tauD']).max():.2e}")
|
|
r_fifo2 = simulate_chain2(3.0, dict(CH2, prio=False))
|
|
print(f" FIFO, theta_A=0.8, theta_B=0.3 (vektet 0.5): Kend={r_fifo2['K'][-1]:.3f} (scen. 2: {res['s2']['K'][-1]:.3f}), "
|
|
f"tau(end)={r_fifo2['tauA'][-1]:.1f}")
|
|
m10 = calibrate_m(dict(CH2, c=10), P)
|
|
m1 = calibrate_m(dict(CH2, c=1), P)
|
|
print(f" kalibrert m: c=1 -> {m1:.4f} (Del 2: {MD:.4f}); c=10 -> {m10:.3f} mnd; "
|
|
f"M/M/10: tau ved ups=0.8: {m10*(1+erlang_L(0.8,10)/(10*0.8)-1+0):.2f}")
|
|
r_c10f = simulate_chain2(3.0, dict(CH2, prio=False, thetaA=0.5, thetaB=0.5, c=10, m=m10))
|
|
print(f" c=10, FIFO: tau_ref={r_c10f['tau_refA']:.2f} mnd")
|
|
|
|
XI_B_UND = 1.0 - 0.2 / (1.0 - CH2["piA"]) # 20 % av alle saker = 33 % av klasse B
|
|
xiB_und = lambda t: XI_B_UND if t >= T_TILTAK else 1.0
|
|
scen2 = {
|
|
"2": ("2 flaskehals, ingen tiltak", dict()),
|
|
"3": ("3 underordne (B: xi_B 1->%.2f)" % XI_B_UND, dict(xiBfun=xiB_und)),
|
|
"4": ("4 utvide", dict(Zfun=Z_utv)),
|
|
"5": ("5 begge", dict(Zfun=Z_utv, xiBfun=xiB_und)),
|
|
"6": ("6 TBT via T (J_D <= 0.95 Z_D)", dict(dbr=0.95)),
|
|
}
|
|
variants = [
|
|
("FIFO, c=1 (referanse, vektet theta)", dict(prio=False)),
|
|
("prioritet, c=1", dict(prio=True)),
|
|
("FIFO, c=10", dict(prio=False, c=10, m=m10)),
|
|
("prioritet, c=10", dict(prio=True, c=10, m=m10)),
|
|
]
|
|
res2 = {}
|
|
for vname, vd in variants:
|
|
print(f" --- {vname} ---")
|
|
for key, (label, sd) in scen2.items():
|
|
ch2 = dict(CH2, **vd, **sd)
|
|
r = simulate_chain2(3.0, ch2)
|
|
res2[(vname, key)] = r
|
|
summary2(label, r)
|
|
r = res2[(vname, "2")]
|
|
print(f" scen. 2: tau_ref A/B = {r['tau_refA']:.2f}/{r['tau_refB']:.2f}; "
|
|
f"tauA@48={r['tauA'][int(48/DT)]:.1f} tauB@48={r['tauB'][int(48/DT)]:.1f}; "
|
|
f"RA_end={r['RA'][-1]:.3f} RB_end={r['RB'][-1]:.3f} frafall B={CH2['varsigma']*r['QB'][-1]:.3f}")
|
|
# figurdata: prioritert kø, scenario 2 og 3
|
|
for key in ("2", "3", "5"):
|
|
rp = res2[("prioritet, c=1", key)]; rf = res2[("FIFO, c=1 (referanse, vektet theta)", key)]
|
|
write_csv(f"data/prio_s{key}.dat", [rp["t"], rp["K"], rp["KA"], rp["KB"], rp["tauA"], rp["tauB"], rf["tauA"], rf["K"], rf["KA"]],
|
|
"t K KA KB tauA tauB tauF KF KAF")
|
|
# alvorlighetsvektet tap: K_A relativt til referanse (0.4) i scen. 2
|
|
print(" tabell (K, K_A/0.4, K_B/0.6, Lambda, tauA, tauB):")
|
|
for vname, _ in variants:
|
|
for key in scen2:
|
|
r = res2[(vname, key)]
|
|
print(f" {vname:36s} scen.{key}: K={r['K'][-1]:.3f} K_A/0.4={r['KA'][-1]/0.4:.3f} K_B/0.6={r['KB'][-1]/0.6:.3f} "
|
|
f"L={r['L'][-1]:.3f} tauA={r['tauA'][-1]:.2f} tauB={r['tauB'][-1]:.2f}")
|
|
rF, rP = res2[("FIFO, c=1 (referanse, vektet theta)", "2")], res2[("prioritet, c=1", "2")]
|
|
K1, L1 = res["s1"]["K"][-1], res["s1"]["L"][-1]
|
|
print(f" merbelastning mot scen. 1 ({K1:.3f}): FIFO {rF['K'][-1]-K1:.3f}, prioritet {rP['K'][-1]-K1:.3f}; "
|
|
f"tillitstap mot scen. 1 ({L1:.3f}): FIFO {L1-rF['L'][-1]:.3f}, prioritet {L1-rP['L'][-1]:.3f}")
|
|
|
|
print("--- (b) regime C (tau=11) med to sakstyper: lånt stabilitet med og uten forrang ---")
|
|
for vname, vd in (("FIFO, c=1", dict(prio=False)), ("prioritet, c=1", dict(prio=True))):
|
|
r = simulate_chain2(11.0, dict(CH2, **vd))
|
|
mx, mn = extrema(r)
|
|
late = r["t"] > 60
|
|
print(f" {vname:16s} scen. 2: amp(K,t>60)={r['K'][late].max()-r['K'][late].min():.3f} topper {mx[:3]} bunner {mn[:3]} "
|
|
f"Lmin={r['L'].min():.3f} tauB(end)={r['tauB'][-1]:.0f}")
|
|
|
|
# --- (c) rekruttering og vippepunkt ---------------------------------------
|
|
print("\n--- (c) rekruttering/smitte: Omega_soc = eps K^p/(K_h^p+K^p), eps=%.2f K_h=%.1f p=%d ---" % (SOC["eps"], SOC["Kh"], SOC["pexp"]))
|
|
Om0s = P["Omega0"] - SOC["eps"] * hill(1.0)
|
|
Kst = SOC["Kh"] * ((SOC["pexp"] - 1) / (SOC["pexp"] + 1)) ** (1 / SOC["pexp"])
|
|
print(f" Omega_0'={Om0s:.4f}; eps*H'(1)={SOC['eps']*hill_d(1.0):.4f} -> a_eff = a - eps H'(1) = {0.12 - SOC['eps']*hill_d(1.0):.4f} (mot 0.12); "
|
|
f"max H' ved K={Kst:.3f}: {hill_d(Kst):.4f}, nødvendig betingelse eps > a/max H' = {0.12/hill_d(Kst):.3f}")
|
|
eq_open = soc_equilibria(P, T=1.0)
|
|
print(f" likevekter, fast T=1 (åpen sløyfe): {np.round(eq_open, 3)} (lav stabil, vippepunkt, høy stabil)")
|
|
for kp_, Tm in ((3.0, 3.0), (1.0, 3.0), (3.0, 1.5), (3.0, 1.2), (3.0, 1.1)):
|
|
eq = soc_equilibria(P, kp=kp_, Tmax=Tm)
|
|
print(f" likevekter med P-regel kp={kp_} og T<= {Tm}: {np.round(eq, 3)}")
|
|
eq = soc_equilibria(dict(P, ksat=1.5), kp=3.0, Tmax=1.5)
|
|
print(f" likevekter med P-regel kp=3, T<=1.5 og politimetning K_sat=1.5: {np.round(eq, 3)}")
|
|
eq = soc_equilibria(dict(P, ksat=1.5), T=1.0)
|
|
print(f" likevekter, fast T=1 og politimetning K_sat=1.5: {np.round(eq, 3)}")
|
|
# tau_c i lav likevekt med rekruttering (a_eff lavere)
|
|
a_soc = 0.12 - SOC["eps"] * hill_d(1.0); b_soc = 0.30
|
|
print(f" lineær analyse i K=1 med rekruttering: a={a_soc:.4f} b={b_soc:.2f} tau_c={critical_delay(a_soc,b_soc):.2f} mnd (mot 7.21)")
|
|
print(" simulering: forbigående sjokk Omega x1.6 i mnd 12-36")
|
|
soc_cases = [
|
|
("fast T=1", "fast", 3.0, P),
|
|
("P kp=3 tau=3", "P", 3.0, P),
|
|
("P kp=3 tau=6", "P", 6.0, P),
|
|
("P kp=3 tau=3, T<=1.5, K_sat=1.5", "P", 3.0, dict(P, Tmax=1.5, ksat=1.5)),
|
|
("P kp=3 tau=6, T<=1.5, K_sat=1.5", "P", 6.0, dict(P, Tmax=1.5, ksat=1.5)),
|
|
]
|
|
cols = []
|
|
for lab, ctrl, tau, pp in soc_cases:
|
|
t, K, U, L, T = simulate_soc(ctrl, tau, pp)
|
|
cols.append(K)
|
|
i36, i48 = int(36 / DT), int(48 / DT)
|
|
print(f" {lab:36s} K@36={K[i36]:.3f} K@48={K[i48]:.3f} Kmax={K.max():.3f} Kend={K[-1]:.3f} "
|
|
f"Tmax={T.max():.2f} Tend={T[-1]:.2f} Lmin={L.min():.3f} Lend={L[-1]:.3f}")
|
|
write_csv("data/vippe.dat", [t] + cols, "t fast p3 p6 p3cap p6cap")
|
|
|
|
# ========================================================================
|
|
print("\n=== Del 4: endogene tiltak -- fra observasjon til beslutning (avsnitt 8.5) ===")
|
|
print("--- grensetilfelle: faste tidsplaner gjenskaper scenario 2 og 5 ---")
|
|
e2 = simulate_endo("sched")
|
|
e5 = simulate_endo("sched", xi_sched=xi_und, Z_sched=Z_utv)
|
|
print(f" sched = scen. 2: max|dK|={np.abs(e2['K'] - res['s2']['K']).max():.2e} "
|
|
f"max|dtauD|={np.abs(e2['tauD'] - res['s2']['tauD']).max():.2e}")
|
|
print(f" sched = scen. 5: max|dK|={np.abs(e5['K'] - res['s5']['K']).max():.2e} "
|
|
f"max|dtauD|={np.abs(e5['tauD'] - res['s5']['tauD']).max():.2e}")
|
|
TEND4 = 180.0
|
|
print(f"--- beslutningsregler, regime A (kp=3, tau=3), {TEND4:.0f} måneder (tabell og figur i 8.5) ---")
|
|
print(f" parametre: lagD={ENDO['lagD']}, xi_min={ENDO['xi_min']}, tau_thr={ENDO['tau_thr']}, k_xi={ENDO['k_xi']:.4f}/mnd, "
|
|
f"tau_xi={ENDO['tau_xi']}, ups*={ENDO['ups_star']}, dødsone [{ENDO['ups_lo']}, {ENDO['ups_hi']}], kZ={ENDO['kZ']}, "
|
|
f"realisering [{ENDO['L1']:.0f}, {ENDO['L2']:.0f}] mnd, tau_down={ENDO['tau_down']}, T_clear={ENDO['T_clear']}, H_pred={ENDO['H_pred']}")
|
|
endo_runs = [
|
|
("e2", "2 flaskehals, ingen tiltak", "sched", dict(), None, None),
|
|
("e5", "5 faste vedtak (mnd 24)", "sched", dict(), xi_und, Z_utv),
|
|
("des", "desentralisert (observert strøm)", "desentral", dict(), None, None),
|
|
("des05", "desentralisert, kZ=0.5", "desentral", dict(kZ=0.5), None, None),
|
|
("desp", "desentralisert med prognose", "desentral_prog", dict(), None, None),
|
|
("ko", "koordinert (T_clear=24)", "koord", dict(), None, None),
|
|
("ko12", "koordinert (T_clear=12)", "koord", dict(T_clear=12.0), None, None),
|
|
("ko48", "koordinert (T_clear=48)", "koord", dict(T_clear=48.0), None, None),
|
|
("koinf", "koordinert (T_clear=inf)", "koord", dict(T_clear=np.inf), None, None),
|
|
("koq", "koordinert + kvalitetssløyfe", "koord", dict(qloop=True), None, None),
|
|
("desq", "desentralisert + kvalitetssløyfe", "desentral", dict(qloop=True), None, None),
|
|
]
|
|
eres = {}
|
|
for key, label, pol, e, xs, zs in endo_runs:
|
|
r = simulate_endo(pol, e=e, tend=TEND4, xi_sched=xs, Z_sched=zs)
|
|
eres[key] = r
|
|
summary_endo(label, r)
|
|
for line in r["log"]:
|
|
print(" " + line)
|
|
for key in ("e5", "des", "desp", "ko"):
|
|
r = eres[key]
|
|
write_csv(f"data/endo_{key}.dat", [r["t"], r["K"], r["L"], r["T"], r["Z"] / ZD0, r["xi"], r["tauD"], r["q"]],
|
|
"t K L T Z xi tauD q")
|
|
# rekkefølge: når reagerer hver sløyfe?
|
|
for key in ("des", "ko"):
|
|
r = eres[key]; t4 = r["t"]
|
|
k_xi = np.argmax(r["xi"] < 0.99); k_xi8 = np.argmax(r["xi"] < 0.81)
|
|
k_z = np.argmax(r["Z"] > ZD0 * 1.001)
|
|
k_rec = np.argmax((t4 > 60) & (r["xi"] > 0.9))
|
|
print(f" {key}: xi<0.99 fra mnd {t4[k_xi]:.0f}, xi<0.81 fra mnd {t4[k_xi8]:.0f}, kapasitet øker fra mnd {t4[k_z]:.0f}, "
|
|
f"xi>0.9 igjen fra mnd {t4[k_rec]:.0f}; J_D ved mnd 33 (budsjettsignal)={r['JD'][int(33 / DT)]:.3f}, "
|
|
f"ups ved mnd 60={r['JD'][int(60 / DT)] / r['Z'][int(60 / DT)]:.2f}")
|
|
r = eres["des"]
|
|
from scipy.signal import argrelextrema
|
|
Zk = r["Z"][::20] / ZD0; tk = r["t"][::20]
|
|
mx = argrelextrema(Zk, np.greater)[0]; mn = argrelextrema(Zk, np.less)[0]
|
|
print(f" des: kapasitetstopper {[(int(tk[i]), round(float(Zk[i]), 2)) for i in mx]} bunner "
|
|
f"{[(int(tk[i]), round(float(Zk[i]), 2)) for i in mn]}; K-topper {extrema(r)[0][:4]}")
|
|
for key in ("des", "ko"):
|
|
r = eres[key]
|
|
print(f" {key}: Z/Z0 og xi ved mnd " + ", ".join(
|
|
f"{m}: {r['Z'][int(m / DT)] / ZD0:.2f}/{r['xi'][int(m / DT)]:.2f}" for m in (36, 48, 60, 72, 96, 108, 120, 138, 150, 180)))
|
|
print(f" ko: tauD topp={eres['ko']['tauD'].max():.1f} ved mnd {eres['ko']['t'][np.argmax(eres['ko']['tauD'])]:.0f}; "
|
|
f"des: tauD topp={eres['des']['tauD'].max():.1f} ved mnd {eres['des']['t'][np.argmax(eres['des']['tauD'])]:.0f}; "
|
|
f"ko: Lend(120)={eres['ko']['L'][int(120 / DT)]:.3f}, koq: Lend(120)={eres['koq']['L'][int(120 / DT)]:.3f}, "
|
|
f"koq: qmax={eres['koq']['q'].max():.3f}")
|
|
print("--- ren kapasitetssløyfe med gitt innstrømning (J=0.134): Z/Z_behov, bestilt kapasitet regnet med / ikke regnet med ---")
|
|
for L2 in (24.0, 36.0):
|
|
for kZ in (0.5, 1.0):
|
|
for sl in (True, False):
|
|
t4, Zr = pure_capacity_loop(kZ, L2=L2, supply_line=sl)
|
|
late = Zr[t4 > 300]; early = Zr[(t4 > 100) & (t4 < 300)]
|
|
print(f" L2={L2:.0f} kZ={kZ}: {'bevilget regnes med ' if sl else 'bare realisert '}: "
|
|
f"Zmax={Zr.max():.2f} amp(100-300)={early.max() - early.min():.3f} amp(>300)={late.max() - late.min():.3f}; "
|
|
f"Hayes a=0, b=kZ/12, tau_Z=(L1+L2)/2: stabil hvis kZ < {np.pi / 2 * 12 / ((6 + L2) / 2):.2f}")
|