function u = ctrl_heatup(t, x, plant, ref)
% CTRL_HEATUP  Ramp T_avg toward a target at a bounded rate.
%
%  Structure:
%       u_ff   = -alpha_f*(T_f - T_f0) - alpha_c*(T_c - T_c0)         % cancel feedback
%       T_ref  = min(ref.T_start + ref.ramp_rate * t, ref.T_target)   % ramped reference
%       u_unsat = u_ff + Kp * (T_ref - T_avg)                         % P on error
%       u      = sat(u_unsat, ref.u_min, ref.u_max)                   % bounded rod worth
%
%  Why saturation:
%    Without it the P gain can push u toward prompt-supercritical as the
%    cold-hot feedback bias unwinds late in the ramp.  Capping u at
%    +0.5*beta guarantees rho_total < beta (below prompt), which in
%    turn bounds the neutron-kinetics excursion rate for reachability.
%
%  IMPORTANT for reach analysis:
%    The sat() is a 3-mode piecewise-affine system (sat-low / linear /
%    sat-high).  Under linear reach assumptions it must either be
%    (a) proven dormant (u_unsat stays in [u_min, u_max] across the
%        reach set — trivial to check, expensive to over-approximate
%        tightly), or
%    (b) handled as an explicit hybrid automaton nested under the DRC
%        mode, with transitions when u_unsat crosses the saturation
%        bounds.
%    The current reach_operation.m assumes (a) implicitly.  Heatup
%    reach will need option (b) because u_unsat is near the +0.5*beta
%    bound during early ramp.
%
%  Why no integrator:
%    Ramp tracking has a structural lag proportional to ramp_rate / Kp_eff.
%    Acceptable because the DRC exits heatup on a predicate window
%    (t_avg_in_range & p_above_crit), not on zero steady-state error.
%    Adding PI would double-count the intrinsic plant integrator
%    (thermal mass) and make anti-windup a hybrid transition.
%
%  IMPORTANT caveat on the cancellation u_ff:
%    The feedback linearization works only if the controller's values
%    of alpha_f, alpha_c match the plant's.  In reality alpha drifts:
%    alpha_f ~20% over burnup, alpha_c by ~10x across boron dilution,
%    plus xenon.  With alpha_true = alpha_nom*(1+delta), the
%    cancellation leaves residual reactivity delta*alpha*dT that the
%    P term cleans up — stabilization still holds, but the clean
%    "rho_total = Kp*e" property is gone.  A robust deployment would
%    treat alpha as an interval and propagate parametric uncertainty
%    through reach (zonotope with parameter generators), or add
%    adaptive alpha estimation.  Idealized here.
%
%  Inputs:
%    t      - time [s]
%    x      - state vector (10 x 1)
%    plant  - parameter struct (alpha_f, alpha_c, T_f0, T_c0 used)
%    ref    - struct with fields:
%               .T_start     starting T_avg [C]
%               .T_target    final T_avg [C]
%               .ramp_rate   desired dT_avg/dt [C/s]
%               .u_min       (optional) lower saturation [dk/k]; default -5*beta
%               .u_max       (optional) upper saturation [dk/k]; default +0.5*beta

    Kp = 1e-4;   % [dk/k per K]

    T_f   = x(8);
    T_c   = x(9);
    T_avg = T_c;

    % Feedforward: cancel intrinsic temperature feedback so rho_total = Kp*e
    % (before saturation).
    u_ff = -plant.alpha_f * (T_f   - plant.T_f0) ...
           -plant.alpha_c * (T_avg - plant.T_c0);

    % Ramped reference, clamped at target.
    T_ref = min(ref.T_start + ref.ramp_rate * t, ref.T_target);

    e        = T_ref - T_avg;
    u_unsat  = u_ff + Kp * e;

    % Saturation bounds (defaults keep rod worth subcritical-prompt).
    if isfield(ref, 'u_min'), u_min = ref.u_min; else, u_min = -5 * plant.beta; end
    if isfield(ref, 'u_max'), u_max = ref.u_max; else, u_max = 0.5 * plant.beta; end

    u = min(max(u_unsat, u_min), u_max);

end