b30e54a8f7
Full-stack app for participatory water pricing using Bezier curves. - Backend: FastAPI + SQLAlchemy + SQLite with JWT auth - Frontend: Nuxt 4 + TypeScript with interactive SVG editor - Math engine: cubic Bezier tarification with Cardano solver - Admin: commune management, household import, vote monitoring, CMS - Citizen: interactive curve editor, vote submission - Docker-compose deployment ready Includes fixes for: - Impact table snake_case/camelCase property mismatch - CMS content backend API + frontend editor (was stub) - Admin route protection middleware - Public content display on commune page - Vote confirmation page link fix Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
119 lines
3.6 KiB
Python
119 lines
3.6 KiB
Python
"""
|
|
Integral computation for Bézier tariff curves.
|
|
|
|
Ported from eau.py:169-211 (NewModel.computeIntegrals).
|
|
Pure Python + numpy, no matplotlib.
|
|
"""
|
|
|
|
import numpy as np
|
|
|
|
|
|
def compute_integrals(
|
|
volume: float,
|
|
vinf: float,
|
|
vmax: float,
|
|
pmax: float,
|
|
a: float,
|
|
b: float,
|
|
c: float,
|
|
d: float,
|
|
e: float,
|
|
) -> tuple[float, float, float]:
|
|
"""
|
|
Compute (alpha1, alpha2, beta2) for a given consumption volume.
|
|
|
|
The total bill for a household consuming `volume` m³ is:
|
|
bill = abo + alpha1 * p0 + alpha2 * p0 + beta2
|
|
|
|
where p0 is the inflection price (computed separately to balance revenue).
|
|
|
|
Args:
|
|
volume: consumption in m³ for this household
|
|
vinf: inflection volume separating the two tiers
|
|
vmax: maximum volume (price = pmax at this volume)
|
|
pmax: maximum price per m³
|
|
a, b: shape parameters for tier 1 Bézier curve
|
|
c, d, e: shape parameters for tier 2 Bézier curve
|
|
|
|
Returns:
|
|
(alpha1, alpha2, beta2) tuple
|
|
"""
|
|
if volume <= vinf:
|
|
# Tier 1 only
|
|
T = _solve_tier1_t(volume, vinf, b)
|
|
alpha1 = _compute_alpha1(T, vinf, a, b)
|
|
return alpha1, 0.0, 0.0
|
|
else:
|
|
# Full tier 1 (T=1) + partial tier 2
|
|
alpha1 = _compute_alpha1(1.0, vinf, a, b)
|
|
|
|
# Tier 2
|
|
wmax = vmax - vinf
|
|
T = _solve_tier2_t(volume - vinf, wmax, c, d)
|
|
|
|
uu = _compute_uu(T, c, d, e)
|
|
alpha2 = (volume - vinf) - 3 * uu * wmax
|
|
beta2 = 3 * pmax * wmax * uu
|
|
return alpha1, alpha2, beta2
|
|
|
|
|
|
def _solve_tier1_t(volume: float, vinf: float, b: float) -> float:
|
|
"""Find T such that v(T) = volume for tier 1."""
|
|
if volume == 0:
|
|
return 0.0
|
|
if volume >= vinf:
|
|
return 1.0
|
|
|
|
# Solve: vinf * [(1 - 3b) * T³ + 3b * T²] = volume
|
|
# => (1-3b) * T³ + 3b * T² - volume/vinf = 0
|
|
p = [1 - 3 * b, 3 * b, 0, -volume / vinf]
|
|
roots = np.roots(p)
|
|
roots = np.unique(roots)
|
|
real_roots = np.real(roots[np.isreal(roots)])
|
|
mask = (real_roots <= 1.0) & (real_roots >= 0.0)
|
|
return float(real_roots[mask][0])
|
|
|
|
|
|
def _solve_tier2_t(w: float, wmax: float, c: float, d: float) -> float:
|
|
"""Find T such that w(T) = w for tier 2, where w = volume - vinf."""
|
|
if w == 0:
|
|
return 0.0
|
|
if w >= wmax:
|
|
return 1.0
|
|
|
|
# Solve: wmax * [(3(c+d-cd)-2)*T³ + 3(1-2c-d+cd)*T² + 3c*T] = w
|
|
p = [
|
|
3 * (c + d - c * d) - 2,
|
|
3 * (1 - 2 * c - d + c * d),
|
|
3 * c,
|
|
-w / wmax,
|
|
]
|
|
roots = np.roots(p)
|
|
roots = np.unique(roots)
|
|
real_roots = np.real(roots[np.isreal(roots)])
|
|
mask = (real_roots <= 1.0 + 1e-10) & (real_roots >= -1e-10)
|
|
if not mask.any():
|
|
# Fallback: closest root to [0,1]
|
|
return float(np.clip(np.real(roots[0]), 0.0, 1.0))
|
|
return float(np.clip(real_roots[mask][0], 0.0, 1.0))
|
|
|
|
|
|
def _compute_alpha1(T: float, vinf: float, a: float, b: float) -> float:
|
|
"""Compute alpha1 coefficient for tier 1."""
|
|
return 3 * vinf * (
|
|
T**6 / 6 * (-9 * a * b + 3 * a + 6 * b - 2)
|
|
+ T**5 / 5 * (24 * a * b - 6 * a - 13 * b + 3)
|
|
+ 3 * T**4 / 4 * (-7 * a * b + a + 2 * b)
|
|
+ T**3 / 3 * 6 * a * b
|
|
)
|
|
|
|
|
|
def _compute_uu(T: float, c: float, d: float, e: float) -> float:
|
|
"""Compute the uu intermediate value for tier 2."""
|
|
return (
|
|
(-3 * c * d + 9 * e * c * d + 3 * c - 9 * e * c + 3 * d - 9 * e * d + 6 * e - 2) * T**6 / 6
|
|
+ (2 * c * d - 15 * e * c * d - 4 * c + 21 * e * c - 2 * d + 15 * e * d - 12 * e + 2) * T**5 / 5
|
|
+ (6 * e * c * d + c - 15 * e * c - 6 * e * d + 6 * e) * T**4 / 4
|
|
+ (3 * e * c) * T**3 / 3
|
|
)
|