gl29/kervent-projections
0
1# optimization.py2"""Financial parameter optimizers for the Kervent simulation.3 4Three modes:5 - compute_min_sci_rent: analytical minimum rent (no scipy, pure arithmetic).6 - optimize_for_breakeven: find parameters that minimize consolidated break-even7 month, using a continuous interpolated metric + scipy differential_evolution.8 - optimize_for_family_income: find parameters that maximize total family treasury9 at a given horizon, same optimizer.10 11All optimizers use copy.deepcopy on simulation inputs to avoid mutating12Streamlit session state.13 14Solvency penalty:15 French cash-flow insolvency test — an entity that cannot pay dues when due is16 insolvent. We penalise every euro-month of deficit across all entities heavily17 so that the optimizer strongly prefers solvent solutions.18"""19 20import copy21from dataclasses import dataclass, field22 23import numpy as np24from scipy.optimize import differential_evolution25 26from core import SimContext, EntitySimulationResult, run_entity_simulation27from ledger import (28 CapExLedger,29 OpExLedger,30 Phase,31 SCIEntity,32 OpCoEntity,33 PersonEntity,34)35 36# ============================================================37# Constants38# ============================================================39 40# Penalty weight per euro-month of insolvency (large enough to dominate the objective)41_SOLVENCY_PENALTY_PER_EURO = 0.142 43 44# ============================================================45# Analytical minimum SCI rent46# ============================================================47 48 49def compute_min_sci_rent(50 phases: list[Phase],51 sci: SCIEntity,52 max_t: int,53 inflation: float,54) -> float:55 """Compute the minimum constant monthly rent (in today's euros) that keeps56 the SCI cash-flow solvent every month.57 58 The rent is configured in today's euros and inflates with the project;59 this function finds the lowest value such that rent * inflation_factor[t]60 covers SCI obligations at every month t.61 62 Only phase-0 (baseline property) loans are considered analytically, since63 phase-2+ trigger months are not known in advance. The result is therefore64 a conservative lower bound — add a safety margin in the UI.65 """66 baseline_capex_arr = sci.capex.compute(max_t)67 baseline_opex_arr = sci.opex.compute(max_t, inflation)68 69 # SCI loan schedule (baseline property only)70 # The downpayment (CapEx out-of-pocket) is an equity injection from savings —71 # NOT a recurring obligation that the rent needs to cover. Only debt service72 # (mortgage) and operating costs are ongoing rent obligations.73 phase_0 = phases[0] if phases else None74 sci_loan_payment = np.zeros(max_t)75 if phase_0 is not None and phase_0.new_loan.loan_years > 0:76 for month, amount in enumerate(baseline_capex_arr):77 if amount <= 0 or phase_0.new_loan.bank_pct <= 0:78 continue79 p, i = phase_0.new_loan.compute_schedule(max_t, amount, start_delay=month)80 sci_loan_payment += p + i81 82 # Total SCI recurring monthly obligations (debt service + operating costs only)83 sci_obligations = baseline_opex_arr + sci_loan_payment84 85 # Rent is in today's euros and inflates: rent_today * infl[t] >= sci_obligations[t]86 # => rent_today >= sci_obligations[t] / infl[t] for all t87 t_arr = np.arange(max_t)88 years_arr = t_arr // 1289 infl_arr = (1.0 + inflation) ** years_arr90 infl_arr = np.where(infl_arr == 0, 1.0, infl_arr)91 92 required_today = sci_obligations / infl_arr93 return float(np.max(required_today))94 95 96# ============================================================97# Simulation wrapper (deep-copy safe)98# ============================================================99 100 101def _run_sim_with_params(102 ctx: SimContext,103 phases: list[Phase],104 base_model_configs: dict,105 sci: SCIEntity,106 opco: OpCoEntity,107 persons: list[PersonEntity],108 params: dict[str, float],109) -> EntitySimulationResult:110 """Run the simulation with parameter overrides.111 112 Uses deep copies of all mutable inputs so the caller's session state is113 never mutated regardless of how many times this is called.114 115 Supported param keys:116 "monthly_rent" -> sci.monthly_rent117 "loan_years_0" -> phases[0].new_loan.loan_years (rounded to int)118 "bank_pct_0" -> phases[0].new_loan.bank_pct119 "target_salary_0" -> persons[0].target_monthly_salary (if persons exist)120 "min_salary_0" -> persons[0].min_monthly_salary121 "min_cash_1" -> phases[1].trigger["min_cash"] (if phase exists)122 "min_dscr_1" -> phases[1].trigger["min_dscr"]123 """124 sci_c = copy.deepcopy(sci)125 opco_c = copy.deepcopy(opco)126 persons_c = copy.deepcopy(persons)127 phases_c = copy.deepcopy(phases)128 129 if "monthly_rent" in params:130 sci_c.monthly_rent = float(params["monthly_rent"])131 if "loan_years_0" in params and phases_c:132 phases_c[0].new_loan.loan_years = max(1, int(round(params["loan_years_0"])))133 if "bank_pct_0" in params and phases_c:134 phases_c[0].new_loan.bank_pct = float(np.clip(params["bank_pct_0"], 0.0, 0.95))135 if "target_salary_0" in params and persons_c:136 persons_c[0].target_monthly_salary = float(params["target_salary_0"])137 if "min_salary_0" in params and persons_c:138 persons_c[0].min_monthly_salary = float(params["min_salary_0"])139 if "min_cash_1" in params and len(phases_c) > 1:140 phases_c[1].trigger["min_cash"] = float(params["min_cash_1"])141 if "min_dscr_1" in params and len(phases_c) > 1:142 phases_c[1].trigger["min_dscr"] = float(params["min_dscr_1"])143 if "initial_capital_sci" in params:144 sci_c.initial_capital = float(params["initial_capital_sci"])145 if "initial_capital_opco" in params:146 opco_c.initial_capital = float(params["initial_capital_opco"])147 148 return run_entity_simulation(149 ctx,150 phases_c,151 base_model_configs,152 sci_c,153 opco_c,154 persons_c,155 )156 157 158# ============================================================159# Solvency penalty + continuous break-even160# ============================================================161 162 163def _solvency_penalty(result: EntitySimulationResult) -> float:164 """Sum of all deficit-euro-months across SCI, OpCo, and persons.165 166 Weighted heavily so the optimizer strongly avoids solutions where any167 entity cannot pay its dues.168 If the simulation failed and produced NaNs, add a massive penalty based on169 how early it failed.170 """171 penalty = 0.0172 for arr in [result.sci_treasury, result.opco_treasury] + list(173 result.person_treasuries.values()174 ):175 nan_mask = np.isnan(arr)176 if np.any(nan_mask):177 t_fail = np.argmax(nan_mask)178 T = len(arr)179 # Massive penalty for dying early180 penalty += (T - t_fail) * 1_000_000.0181 valid_arr = arr[:t_fail]182 penalty += float(np.sum(np.abs(np.minimum(valid_arr, 0.0)))) * _SOLVENCY_PENALTY_PER_EURO183 else:184 penalty += float(np.sum(np.abs(np.minimum(arr, 0.0)))) * _SOLVENCY_PENALTY_PER_EURO185 return penalty186 187 188def _continuous_breakeven(result: EntitySimulationResult) -> float:189 """Continuous break-even metric via linear interpolation of the zero-crossing.190 191 Avoids the discrete step-function problem that defeats gradient-based192 optimizers when using the integer break_even_idx directly.193 194 Returns a float in [0, max_t + large_penalty]. Smaller is better.195 """196 consol = result.consolidated_treasury197 max_t = len(consol)198 199 if result.failure_month is not None:200 # Massive penalty if it didn't even survive201 return float(max_t) + (max_t - result.failure_month) * 1000.0202 203 trough_idx = int(np.nanargmin(consol))204 post_trough = consol[trough_idx:]205 206 if not np.any(post_trough >= 0):207 # Never recovers: return max_t plus a penalty proportional to final deficit208 return float(max_t) + float(np.abs(consol[-1]))209 210 cross_rel = int(np.argmax(post_trough >= 0))211 be_idx = trough_idx + cross_rel212 213 if be_idx == 0:214 return 0.0215 216 val_before = consol[be_idx - 1]217 val_after = consol[be_idx]218 if val_after != val_before:219 fraction = -val_before / (val_after - val_before)220 else:221 fraction = 0.0222 223 return float(be_idx - 1) + float(np.clip(fraction, 0.0, 1.0))224 225 226# ============================================================227# Optimization result228# ============================================================229 230 231@dataclass232class OptimizationResult:233 objective: str234 optimal_params: dict[str, float]235 optimal_value: float236 baseline_value: float237 converged: bool238 message: str239 n_evaluations: int240 # Full simulation result at optimal parameters (for before/after charts)241 optimal_simulation: EntitySimulationResult | None = None242 243 244# ============================================================245# Shared optimizer core246# ============================================================247 248# Human-readable labels for free variable keys249PARAM_LABELS: dict[str, str] = {250 "monthly_rent": "Loyer mensuel SCI (€/mois)",251 "loan_years_0": "Durée du prêt Phase 1 (ans)",252 "bank_pct_0": "Part bancaire Phase 1 (%)",253 "target_salary_0": "Salaire cible famille (€/mois)",254 "min_salary_0": "Salaire garanti famille (€/mois)",255 "min_cash_1": "Trésorerie min. déclenchement Phase 2 (€)",256 "min_dscr_1": "DSCR min. déclenchement Phase 2",257 "initial_capital_sci": "Apport Initial SCI (€)",258 "initial_capital_opco": "Apport Initial OpCo (€)",259}260 261# Default bounds for each parameter key262DEFAULT_BOUNDS: dict[str, tuple[float, float]] = {263 "monthly_rent": (0.0, 8000.0),264 "loan_years_0": (5.0, 30.0),265 "bank_pct_0": (0.0, 0.95),266 "target_salary_0": (0.0, 5000.0),267 "min_salary_0": (0.0, 3000.0),268 "min_cash_1": (0.0, 50000.0),269 "min_dscr_1": (0.5, 3.0),270 "initial_capital_sci": (0.0, 200_000.0),271 "initial_capital_opco": (0.0, 100_000.0),272}273 274 275def _build_objective(276 objective_type: str,277 ctx: SimContext,278 phases: list[Phase],279 base_model_configs: dict,280 sci: SCIEntity,281 opco: OpCoEntity,282 persons: list[PersonEntity],283 free_var_keys: list[str],284 horizon_months: int | None = None,285):286 """Return a scalar objective function f(x) -> float for scipy.287 288 x is a 1-D array of free variable values in the order of free_var_keys.289 """290 call_count = [0]291 292 def objective(x: np.ndarray) -> float:293 call_count[0] += 1294 params = dict(zip(free_var_keys, x))295 result = _run_sim_with_params(296 ctx,297 phases,298 base_model_configs,299 sci,300 opco,301 persons,302 params,303 )304 305 penalty = _solvency_penalty(result)306 307 if objective_type == "min_breakeven":308 primary = _continuous_breakeven(result)309 elif objective_type == "max_family_income":310 # Maximise total family treasury at end of horizon311 h = horizon_months or ctx.max_t312 h = min(h, ctx.max_t)313 family_total = sum(arr[h - 1] for arr in result.person_treasuries.values())314 primary = -family_total # negate because we minimise315 elif objective_type == "min_capital":316 # Minimize the required initial capital (assuming the params map correctly)317 primary = float(params.get("initial_capital_sci", getattr(sci, "initial_capital", 0.0))) + \318 float(params.get("initial_capital_opco", getattr(opco, "initial_capital", 0.0)))319 else:320 raise ValueError(f"Unknown objective: {objective_type}")321 322 return primary + penalty323 324 return objective, call_count325 326 327def _run_optimization(328 objective_type: str,329 ctx: SimContext,330 phases: list[Phase],331 base_model_configs: dict,332 sci: SCIEntity,333 opco: OpCoEntity,334 persons: list[PersonEntity],335 free_var_keys: list[str],336 bounds_override: dict[str, tuple[float, float]] | None = None,337 horizon_months: int | None = None,338 maxiter: int = 200,339 popsize: int = 12,340 seed: int = 42,341) -> OptimizationResult:342 """Run differential_evolution and return an OptimizationResult."""343 344 bo = bounds_override or {}345 bounds = [bo.get(k, DEFAULT_BOUNDS.get(k, (0.0, 1.0))) for k in free_var_keys]346 347 obj_fn, call_count = _build_objective(348 objective_type,349 ctx,350 phases,351 base_model_configs,352 sci,353 opco,354 persons,355 free_var_keys,356 horizon_months,357 )358 359 # Baseline value (current parameters)360 baseline_params = {k: 0.0 for k in free_var_keys} # will read from current inputs361 baseline_result = _run_sim_with_params(362 ctx,363 phases,364 base_model_configs,365 sci,366 opco,367 persons,368 baseline_params,369 )370 if objective_type == "min_breakeven":371 baseline_value = _continuous_breakeven(baseline_result)372 elif objective_type == "min_capital":373 baseline_value = float(getattr(sci, "initial_capital", 0.0) + getattr(opco, "initial_capital", 0.0))374 else:375 h = horizon_months or ctx.max_t376 h = min(h, ctx.max_t)377 baseline_value = sum(378 arr[h - 1] for arr in baseline_result.person_treasuries.values()379 )380 381 de_result = differential_evolution(382 obj_fn,383 bounds=bounds,384 maxiter=maxiter,385 popsize=popsize,386 seed=seed,387 tol=1e-4,388 workers=1, # Streamlit is not fork-safe; keep single-threaded389 polish=True,390 )391 392 optimal_params = dict(zip(free_var_keys, de_result.x))393 394 # Run one final clean simulation at optimal params for charts395 optimal_sim = _run_sim_with_params(396 ctx,397 phases,398 base_model_configs,399 sci,400 opco,401 persons,402 optimal_params,403 )404 405 # Report objective without penalty for display406 if objective_type == "min_breakeven":407 optimal_value = _continuous_breakeven(optimal_sim)408 elif objective_type == "min_capital":409 optimal_value = optimal_params.get("initial_capital_sci", getattr(sci, "initial_capital", 0.0)) + \410 optimal_params.get("initial_capital_opco", getattr(opco, "initial_capital", 0.0))411 else:412 h = horizon_months or ctx.max_t413 h = min(h, ctx.max_t)414 optimal_value = sum(415 arr[h - 1] for arr in optimal_sim.person_treasuries.values()416 )417 418 return OptimizationResult(419 objective=objective_type,420 optimal_params=optimal_params,421 optimal_value=float(optimal_value),422 baseline_value=float(baseline_value),423 converged=bool(de_result.success),424 message=de_result.message,425 n_evaluations=call_count[0],426 optimal_simulation=optimal_sim,427 )428 429 430# ============================================================431# Public API432# ============================================================433 434 435def optimize_for_breakeven(436 ctx: SimContext,437 phases: list[Phase],438 base_model_configs: dict,439 sci: SCIEntity,440 opco: OpCoEntity,441 persons: list[PersonEntity],442 free_var_keys: list[str],443 bounds_override: dict[str, tuple[float, float]] | None = None,444 maxiter: int = 200,445 popsize: int = 12,446) -> OptimizationResult:447 """Find parameter values that minimise consolidated break-even month.448 449 Uses a continuous interpolated break-even proxy and scipy450 differential_evolution, which handles bounds and discontinuous objectives451 better than gradient-based methods.452 """453 return _run_optimization(454 "min_breakeven",455 ctx,456 phases,457 base_model_configs,458 sci,459 opco,460 persons,461 free_var_keys=free_var_keys,462 bounds_override=bounds_override,463 maxiter=maxiter,464 popsize=popsize,465 )466 467 468def optimize_for_family_income(469 ctx: SimContext,470 phases: list[Phase],471 base_model_configs: dict,472 sci: SCIEntity,473 opco: OpCoEntity,474 persons: list[PersonEntity],475 free_var_keys: list[str],476 bounds_override: dict[str, tuple[float, float]] | None = None,477 horizon_months: int | None = None,478 maxiter: int = 200,479 popsize: int = 12,480) -> OptimizationResult:481 """Find parameter values that maximise total family treasury at the horizon.482 483 Maximising means minimising the negative family treasury. Solvency484 violations are penalised so the optimizer prefers solvent solutions.485 """486 return _run_optimization(487 "max_family_income",488 ctx,489 phases,490 base_model_configs,491 sci,492 opco,493 persons,494 free_var_keys=free_var_keys,495 bounds_override=bounds_override,496 horizon_months=horizon_months,497 maxiter=maxiter,498 popsize=popsize,499 )500 501 502def optimize_for_min_capital(503 ctx: SimContext,504 phases: list[Phase],505 base_model_configs: dict,506 sci: SCIEntity,507 opco: OpCoEntity,508 persons: list[PersonEntity],509 free_var_keys: list[str],510 bounds_override: dict[str, tuple[float, float]] | None = None,511 maxiter: int = 200,512 popsize: int = 12,513) -> OptimizationResult:514 """Find parameter values that minimize the required initial capital for survival.515 """516 return _run_optimization(517 "min_capital",518 ctx,519 phases,520 base_model_configs,521 sci,522 opco,523 persons,524 free_var_keys=free_var_keys,525 bounds_override=bounds_override,526 maxiter=maxiter,527 popsize=popsize,528 )529 