ESC Dynamics with Thermal Derating (8-state)

def _duty_command(t: float) -> float:
    """Simulated flight controller PWM commands."""
    if t < 0.5:
        return 0.3          # idle / spin-up
    elif t < 2.0:
        return 0.65         # hover throttle
    elif t < 4.0:
        return 0.85         # climb / aggressive manoeuvre
    elif t < 6.0:
        return 0.65         # back to hover
    elif t < 8.0:
        return 0.4          # descent
    else:
        return 0.1          # motor brake / idle

def _motor_omega(duty: float) -> float:
    """Approximate motor speed from duty cycle for back-EMF."""
    return duty * _V_BUS / _KE * 0.85

def rhs_esc_dynamics(t, y):
    """8-state ESC dynamics with current limiting and thermal derating."""
    dy = np.zeros(8)

    V_out = y[0]
    I_out = y[1]
    T_mos = y[2]
    T_brd = y[3]
    duty = y[4]
    I_cap = y[5]
    V_cap = y[6]

    # PWM filter: track commanded duty with first-order lag
    d_cmd = _duty_command(t)

    # Thermal derating
    if T_mos > _T_DERATE:
        derate = max(0.0, 1.0 - (T_mos - _T_DERATE) / (_T_SHUTDOWN - _T_DERATE))
    else:
        derate = 1.0
    d_limited = d_cmd * derate

    # Current limiting (soft)
    if I_out > _I_SOFT:
        i_limit_factor = max(0.1, 1.0 - (I_out - _I_SOFT) / (_I_MAX - _I_SOFT))
        d_limited *= i_limit_factor

    dy[4] = (d_limited - duty) / _TAU_PWM

    # Output voltage
    omega = _motor_omega(duty)
    V_bemf = _KE * omega
    V_applied = V_cap * duty

    dy[0] = (V_applied - V_out) / (_TAU_PWM * 10)

    # Output current: V = I*R + L*dI/dt + V_bemf
    if _L_MOTOR > 0:
        dy[1] = (V_out - I_out * _R_MOTOR - V_bemf) / _L_MOTOR
    else:
        dy[1] = 0.0

    # Clamp current
    if I_out >= _I_MAX and dy[1] > 0:
        dy[1] = 0.0

    # MOSFET conduction + switching losses
    P_cond = I_out**2 * _R_DS_ON * 2  # 2 FETs in current path
    P_sw = 0.5 * V_cap * I_out * 20e-9 * _F_SW * 2  # switching loss estimate
    P_total = P_cond + P_sw

    # Thermal dynamics: junction
    dy[2] = (P_total - (T_mos - T_brd) / _R_TH_JC) / _C_TH_J

    # Thermal dynamics: board
    dy[3] = ((T_mos - T_brd) / _R_TH_JC - (T_brd - _T_AMB) / _R_TH_BA) / _C_TH_B

    # DC bus capacitor
    I_bus = I_out * duty
    dy[5] = (I_bus - I_cap) / (_TAU_PWM * 5)
    dy[6] = (I_cap - I_bus) / _C_BUS if _C_BUS > 0 else 0.0

    # Cumulative energy dissipated
    dy[7] = P_total

    return dy