"""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. (9).""" 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. (9) 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. (25). # (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 2, 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}")