Kjell Tore Guttormsen
4a9170c3b7
Four references still used the numbering from before the note was extended: the crime dynamics are equation (12), the court queue is equation (28) and the parameter table is table 3. Comments and one printed heading only. Verified by running sim.py before and after: 22/22 data files byte-identical, output differs in that one heading. 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}")
|