feat: add charge_profile for SoC-dependent charge power limits

openapi.yaml

@@ -269,6 +269,39 @@ components:
           maximum: 2
           default: 0
           description: Charging and discharging priority 0..2 compared to other batteries. 2 = highest priority.
+        charge_profile:
+          type: string
+          description: |
+            Named charge profile for SoC-dependent charge power limiting
+            (CC-CV taper). When set, the MILP solver uses iterative constraint
+            tightening to enforce realistic charge power limits at high SoC.
+            Available profiles: lifepo4_conservative, lifepo4_moderate.
+            Individual parameters can be overridden with charge_knee,
+            charge_k, charge_c_rate_float.
+          example: lifepo4_conservative
+        charge_knee:
+          type: number
+          minimum: 0
+          maximum: 100
+          description: |
+            SoC percentage where charge power tapering begins.
+            With charge_profile: overrides the profile default.
+            Without charge_profile: requires charge_k to build a custom profile.
+          example: 65
+        charge_k:
+          type: number
+          minimum: 0
+          description: |
+            Exponential decay constant for the CC-CV taper curve.
+            Higher values mean steeper power reduction above the knee.
+          example: 0.052
+        charge_c_rate_float:
+          type: number
+          minimum: 0
+          description: |
+            Minimum C-rate at the tail end of charging (float/balancing phase).
+            Default: 0.01 (C/100).
+          example: 0.01
 
     TimeSeries:
       type: object

src/optimizer/app.py

@@ -5,6 +5,7 @@
 from flask_restx import Api, Resource, fields
 from werkzeug.exceptions import BadRequest
 
+from .charging_profiles import get_profile
 from .optimizer import BatteryConfig, GridConfig, OptimizationStrategy, Optimizer, TimeSeriesData
 
 app = Flask(__name__)
@@ -80,7 +81,11 @@ def handle_validation_error(error):
     'c_max': fields.Float(required=True, description='Maximum charge power (W)'),
     'd_max': fields.Float(required=True, description='Maximum discharge power (W)'),
     'p_a': fields.Float(required=True, description='Monetary value per Wh at end of the optimization horizon'),
-    'c_priority': fields.Integer(required=False, description='Charging and discharging priority compared to other batteries. 2 = highest priority.')
+    'c_priority': fields.Integer(required=False, description='Charging and discharging priority compared to other batteries. 2 = highest priority.'),
+    'charge_profile': fields.String(required=False, description='Named CC-CV taper profile (e.g. lifepo4_conservative).'),
+    'charge_knee': fields.Float(required=False, description='SoC% where taper begins. Overrides profile default.'),
+    'charge_k': fields.Float(required=False, description='Exponential decay constant for the CC-CV taper curve.'),
+    'charge_c_rate_float': fields.Float(required=False, description='Minimum C-rate at float / tail end of charging (default 0.01).'),
 })
 
 time_series_model = api.model('TimeSeries', {
@@ -157,20 +162,50 @@ def post(self):
             # Parse battery configurations
             batteries = []
             for bat_data in data['batteries']:
+                # Build optional CC-CV taper profile
+                profile = None
+                profile_name = bat_data.get('charge_profile')
+                charge_knee = bat_data.get('charge_knee')
+                charge_k = bat_data.get('charge_k')
+                charge_c_rate_float = bat_data.get('charge_c_rate_float')
+
+                s_capacity = bat_data.get('s_capacity', bat_data['s_max'])
+                c_max = bat_data['c_max']
+
+                if profile_name:
+                    # Named profile with optional overrides
+                    overrides = {}
+                    if charge_knee is not None:
+                        overrides['knee'] = charge_knee
+                    if charge_k is not None:
+                        overrides['k'] = charge_k
+                    if charge_c_rate_float is not None:
+                        overrides['c_rate_float'] = charge_c_rate_float
+                    profile = get_profile(profile_name, **overrides)
+                elif charge_knee is not None and charge_k is not None:
+                    # Custom profile: derive c_rate_max from c_max / s_capacity
+                    profile = {
+                        'c_rate_max': c_max / s_capacity if s_capacity > 0 else 0.25,
+                        'knee': charge_knee,
+                        'k': charge_k,
+                        'c_rate_float': charge_c_rate_float if charge_c_rate_float is not None else 0.01,
+                    }
+
                 batteries.append(BatteryConfig(
                     charge_from_grid=bat_data.get('charge_from_grid', False),
                     discharge_to_grid=bat_data.get('discharge_to_grid', False),
-                    s_capacity=bat_data.get('s_capacity', bat_data['s_max']),
+                    s_capacity=s_capacity,
                     s_min=bat_data['s_min'],
                     s_max=bat_data['s_max'],
                     s_initial=bat_data['s_initial'],
                     p_demand=bat_data.get('p_demand'),
                     s_goal=bat_data.get('s_goal'),
                     c_min=bat_data['c_min'],
-                    c_max=bat_data['c_max'],
+                    c_max=c_max,
                     d_max=bat_data['d_max'],
                     p_a=bat_data['p_a'],
                     c_priority=bat_data.get('c_priority', 0),
+                    charge_profile=profile,
                 ))
 
             # Parse time series data

src/optimizer/charging_profiles.py

@@ -0,0 +1,132 @@
+"""
+Battery charging profiles by cell chemistry.
+
+Each profile defines how the maximum charge power decreases as a function
+of State of Charge (SoC). This models the real-world CC-CV charging behavior
+where the Battery Management System (BMS) reduces the charge current at
+higher SoC to protect cells and minimize resistive losses.
+
+Usage:
+    from charging_profiles import get_profile, max_charge_power
+
+    profile = get_profile("lifepo4_conservative")
+    watts = max_charge_power(profile, capacity_wh=31000, soc_percent=75)
+
+To find values for your own system:
+    - Set c_rate_max to the highest C-rate your battery sustains
+      without the BMS reducing the charge current
+    - Set knee to the SoC percentage where the BMS starts
+      reducing the charge current below c_rate_max
+    - Leave k and c_rate_float at their defaults unless you have
+      your own measurements (fit with Wolfram Alpha:
+      "exponential fit {SoC1,I1}, {SoC2,I2}, ...")
+
+General advice: always stay slightly below your measured values.
+If your BMS starts tapering at 66%, set knee to 65%. If your
+inverter sustains 160A, set c_rate_max to match C/4 or lower.
+This avoids triggering BMS intervention, reduces resistive losses,
+and extends cell life.
+"""
+
+import math
+
+# ---------------------------------------------------------------------------
+# Chemistry profiles
+# ---------------------------------------------------------------------------
+
+CHEMISTRY_PROFILES = {
+    "lifepo4_conservative": {
+        # Charging profile for LiFePO4 cells optimized for longevity and
+        # minimal resistive losses. Derived from real-world measurements
+        # on a mixed environment of aged LiYFePO4 cells (200Ah + 90Ah)
+        # and newer LiFePO4 cells (314Ah) with Victron inverters.
+        #
+        # The charging curve has two phases:
+        #   1. Below 'knee': constant power at 'c_rate_max' rate
+        #   2. Above 'knee': power drops exponentially toward 'c_rate_float'
+        #
+        # Formula above knee:
+        #   P = max(c_rate_float, c_rate_max * e^(-k * (SoC - knee)))
+
+        "c_rate_max": 0.25,
+        # Maximum charge rate below knee, as fraction of capacity.
+        # C/4 means a 31 kWh battery charges at up to 7750 W.
+
+        "knee": 65,
+        # SoC percentage where charge current begins to taper.
+
+        "k": 0.052,
+        # Exponential decay constant (R^2 = 0.996 from BMS measurements).
+
+        "c_rate_float": 0.01,
+        # Minimum charge rate (maintenance/balancing phase).
+        # C/100 means a 31 kWh battery trickle-charges at 310 W.
+    },
+
+    "lifepo4_moderate": {
+        # Moderate profile for newer LiFePO4 cells in good condition.
+
+        "c_rate_max": 0.33,     # C/3
+        "knee": 75,
+        "k": 0.045,
+        "c_rate_float": 0.02,   # C/50
+    },
+}
+
+
+# ---------------------------------------------------------------------------
+# Public API
+# ---------------------------------------------------------------------------
+
+def list_profiles():
+    """Return list of available profile names."""
+    return list(CHEMISTRY_PROFILES.keys())
+
+
+def get_profile(name, **overrides):
+    """
+    Get a charging profile by name with optional overrides.
+
+    Args:
+        name: Profile name (e.g. "lifepo4_conservative")
+        **overrides: Override individual values (e.g. knee=70, c_rate_max=0.2)
+
+    Returns:
+        dict with keys: c_rate_max, knee, k, c_rate_float
+    """
+    if name not in CHEMISTRY_PROFILES:
+        available = ", ".join(CHEMISTRY_PROFILES.keys())
+        raise ValueError(f"Unknown profile '{name}'. Available: {available}")
+
+    profile = dict(CHEMISTRY_PROFILES[name])
+    for key, value in overrides.items():
+        if key not in profile:
+            raise ValueError(f"Unknown parameter '{key}'. Valid: {', '.join(profile.keys())}")
+        profile[key] = value
+
+    return profile
+
+
+def max_charge_power(profile, capacity_wh, soc_percent):
+    """
+    Calculate maximum charge power in watts for a given SoC.
+
+    Args:
+        profile: Charging profile dict from get_profile()
+        capacity_wh: Battery capacity in Wh
+        soc_percent: Current state of charge (0-100)
+
+    Returns:
+        Maximum charge power in watts
+    """
+    c_rate_max = profile["c_rate_max"]
+    knee = profile["knee"]
+    k = profile["k"]
+    c_rate_float = profile["c_rate_float"]
+
+    if soc_percent <= knee:
+        c_rate = c_rate_max
+    else:
+        c_rate = max(c_rate_float, c_rate_max * math.exp(-k * (soc_percent - knee)))
+
+    return capacity_wh * c_rate

src/optimizer/optimizer.py

@@ -1,10 +1,12 @@
+import math
 from dataclasses import dataclass
 from tempfile import TemporaryDirectory
 from typing import Dict, List, Optional
 
 import numpy as np
 import pulp
 
+from .charging_profiles import max_charge_power
 from .settings import OptimizerSettings
 
 
@@ -36,6 +38,7 @@ class BatteryConfig:
     p_demand: Optional[List[float]] = None  # Minimum charge demand (Wh)
     s_goal: Optional[List[float]] = None  # Goal state of charge (Wh)
     c_priority: int = 0
+    charge_profile: Optional[dict] = None  # CC-CV taper profile (charging_profiles format)
 
 
 @dataclass
@@ -115,6 +118,7 @@ def create_model(self):
         self._setup_target_function()
         self._add_energy_balance_constraints()
         self._add_battery_constraints()
+        self._add_charge_profile_constraints()
 
     def _setup_variables(self):
         """
@@ -491,10 +495,80 @@ def _add_battery_constraints(self):
                 # Charge constraint
                 self.problem += self.variables['c'][i][t] <= self.M * (1 - self.variables['z_cd'][i][t])
 
+    def _add_charge_profile_constraints(self):
+        """Add piecewise-linear upper bounds on charge power from CC-CV taper profiles.
+
+        For each battery with a charge_profile, the taper curve P_max(SoC) is
+        a convex decreasing function above the knee. We approximate it with
+        tangent lines at 4 breakpoints. Each tangent is a valid linear upper
+        bound (convexity guarantees the tangent never exceeds the curve), so
+        the MILP enforces realistic charge power limits in a single solve.
+
+        The constraint for slot t uses s[i][t-1] (SoC at slot start). For
+        t=0, s_initial is a constant so the bound is applied directly to
+        the variable's upper bound.
+        """
+        for i, bat in enumerate(self.batteries):
+            if bat.charge_profile is None:
+                continue
+
+            profile = bat.charge_profile
+            cap = bat.s_capacity
+            if cap <= 0:
+                continue
+
+            knee = profile['knee']
+            k = profile['k']
+            c_rate_max = profile['c_rate_max']
+
+            # 4 breakpoints from knee to 100% SoC
+            n_points = 4
+            breakpoints = []
+            for j in range(n_points):
+                soc_pct = knee + (100.0 - knee) * j / (n_points - 1)
+                s_wh = soc_pct / 100.0 * cap
+                p_w = max_charge_power(profile, cap, soc_pct)
+                # derivative of P_max w.r.t. SoC in Wh:
+                # dP/ds = dP/d(soc%) * d(soc%)/ds
+                #       = (-k * c_rate_max * cap * e^(-k*(soc%-knee))) * (100/cap)
+                #       = -k * c_rate_max * 100 * e^(-k*(soc%-knee))
+                dp_ds = -k * c_rate_max * 100.0 * math.exp(-k * (soc_pct - knee))
+                breakpoints.append((s_wh, p_w, dp_ds))
+
+            for s_b, p_b, dp_ds in breakpoints:
+                for t in self.time_steps:
+                    dt_h = self.time_series.dt[t] / 3600.0
+
+                    if t == 0:
+                        # s_initial is a constant: compute a fixed upper bound
+                        p_limit = p_b + dp_ds * (bat.s_initial - s_b)
+                        e_limit = p_limit * dt_h
+                        e_var_ub = bat.c_max * dt_h
+                        # only add if tighter than the existing variable bound
+                        if e_limit < e_var_ub - 1.0:
+                            self.problem += (
+                                self.variables['c'][i][t] <= e_limit,
+                                f"taper_{i}_{t}_bp{breakpoints.index((s_b, p_b, dp_ds))}",
+                            )
+                    else:
+                        # linear constraint: c[i][t] <= (p_b - dp_ds*s_b)*dt_h + dp_ds*dt_h * s[i][t-1]
+                        s_prev = self.variables['s'][i][t - 1]
+                        intercept = (p_b - dp_ds * s_b) * dt_h
+                        slope = dp_ds * dt_h
+                        self.problem += (
+                            self.variables['c'][i][t] <= intercept + slope * s_prev,
+                            f"taper_{i}_{t}_bp{breakpoints.index((s_b, p_b, dp_ds))}",
+                        )
+
     def solve(self) -> Dict:
         """
         Creates the MILP model if none exists and solves the optimization problem.
-        Returns a dictionary with the optimization results
+
+        When batteries have a charge_profile (CC-CV taper), piecewise-linear
+        upper bounds on charge power are part of the model, so a single solve
+        enforces realistic SoC-dependent charge limits.
+
+        Returns a dictionary with the optimization results.
         """
 
         if self.problem is None: