% =========================================================================
% LI AUTO L9: INTEGRATED APU ENERGY STORAGE OPTIMIZATION (EMPIRICAL BATTERY)
% OBJECTIVE: Maximize Total WLTC Range
% CONSTRAINTS: 230km EV Range, 120km/h @ 6% Incline, 3.6 m^2 Powertrain Area
% FIXED PARAMETERS: APU System = 115 kW, Motor = 250 kW
% =========================================================================
clear; clc; close all;

%% 1. PHYSICAL CONSTANTS & ENVIRONMENTAL FACTORS
data.g = 9.81; data.f = 0.01; data.Cd = 0.30; data.A_front = 2.5; 
data.eta_t = 0.94; data.eta_eng = 0.445; data.eta_gen = 0.95; 
data.rho_Ex = 9.5;          % Fuel Energy Density (kWh/L)
data.DoD = 0.93;            % Usable Battery Depth of Discharge
data.P_aux = 3.0;           % Auxiliary system load (kW)

%% 2. FIXED POWERTRAIN SPECS & MERGED DIMENSIONS
data.P_eg_fixed = 115;      % Fixed Merged Engine-Generator System Power (kW)
data.P_motor_fixed = 250;   % Fixed Traction Motor Power (kW)

% Component Weight Scaling (kg per unit)
data.g_f = 1.0;          % kg/L (Fuel + Tank)
data.g_b_slope = 6.16;   % kg/kWh (Battery cell density regression)
data.g_b_fixed = 83.87;  % kg (Fixed pack casing & BMS regression)
data.g_eg = 2.0;         % kg/kW (Combined Engine + Generator System)
data.g_m = 0.5;          % kg/kW (Traction Motor)

% Factory Baseline Reference 
data.m_base = 2755;         
data.V_base = 65;           
data.E_base = 61.0;         % NEW REFERENCE: 61 kWh Battery
data.P_eng_base = 115;      
data.P_motor_base = 420;    

%% 3. CONSOLIDATED 2D COMPONENT AREAS (m^2)
data.A_tank_base = 0.422;   % 1000 * 422 mm
data.A_batt_base = 2.59;    % NEW REFERENCE: Battery Area for 61 kWh

% --- DYNAMIC FRONT/REAR MOTOR AREA SCALING ---
P_rear_base  = 275;         % kW (Rear Motor Base Power)
A_rear_base  = 0.261;       % m^2 (Rear Motor Base Area)
P_front_base = 145;         % kW (Front Motor Base Power)

% Calculate Front Motor Area using the Rear Motor's area-to-power ratio
A_front_base = A_rear_base * (P_front_base / P_rear_base); 

% Total Base Motor Area (Front + Rear representing the full 420 kW)
data.A_motor_base = A_rear_base + A_front_base; 

% Engine + Generator Area (Total 1.0 m^2 front area minus the front motor)
data.A_eg_base = 1.0 - A_front_base; 

% Updated Optimization Targets
data.A_max_limit = 3.550;    % Expanded to 3.60 m^2 to prevent physical solver failure
data.min_ev_range = 200;    % Target EV Range (km)

%% 4. EXTRACT WEIGHT-INDEPENDENT GLIDER MASS
% Calculate baseline battery mass using regression
m_batt_base = (data.g_b_slope * data.E_base) + data.g_b_fixed;

data.m_glider = data.m_base - (data.g_f*data.V_base + m_batt_base + ...
                data.g_eg*data.P_eng_base + data.g_m*data.P_motor_base);

% Calibration scalar for pure electric consumption tracking
data.v_scalar = 1.1052; 

%% 5. EXECUTE NON-LINEAR OPTIMIZATION (fmincon)
% Setup Optimization Variables: x = [Fuel(L), Battery(kWh)]
x0 = [100, 70]; % Initial guess

% Engineering Bounds [LB, UB]
lb = [10,  20];  % Min: 10 Liters, 20 kWh
ub = [200, 180]; % Max: 200 Liters, 180 kWh

% Optimization Options
options = optimoptions('fmincon', 'Display', 'final', 'Algorithm', 'sqp', 'MaxFunctionEvaluations', 2000);

% Run Solver: minimize negative range -> maximizes positive range
[x_opt, neg_max_range] = fmincon(@(x) obj_range(x, data), x0, [], [], [], [], lb, ub, @(x) dynamic_constraints(x, data), options);

max_range = -neg_max_range;

%% 6. POST-PROCESSING FOR PARAMETER PRINTOUT
V_fuel_opt = x_opt(1); E_batt_opt = x_opt(2);

% Final Component Masses
mass_fuel  = data.g_f * V_fuel_opt;
mass_batt  = (data.g_b_slope * E_batt_opt) + data.g_b_fixed;
mass_eg    = data.g_eg * data.P_eg_fixed;
mass_motor = data.g_m * data.P_motor_fixed;
opt_mass   = data.m_glider + mass_fuel + mass_batt + mass_eg + mass_motor;

% Final Component Areas (Using the New Empirical Battery Equation)
A_tank_opt  = data.A_tank_base * (V_fuel_opt / data.V_base);
A_batt_opt  = (E_batt_opt + 1.42) / 23.30;   % <--- EMPIRICAL EQUATION APPLIED HERE
A_eg_opt    = data.A_eg_base   * (data.P_eg_fixed / data.P_eng_base);
A_motor_opt = data.A_motor_base * (data.P_motor_fixed / data.P_motor_base);
opt_Area    = A_tank_opt + A_batt_opt + A_eg_opt + A_motor_opt;

% Final Consumption Dynamics
v_norm = 100;
P_roll_opt = (opt_mass * data.g * data.f * v_norm) / 3600;
P_req_opt = (P_roll_opt + (data.Cd * data.A_front * v_norm^3) / 76140) / data.eta_t;
Ec_opt = P_req_opt * 1; 
Ec_EV_opt = Ec_opt + data.P_aux;
R_EV_opt = ((E_batt_opt * data.DoD) / Ec_EV_opt) * 100;

%% 7. COMPREHENSIVE PARAMETER PRINTOUT
fprintf('\n====================================\n');
fprintf('       POWERTRAIN SPECIFICATIONS        \n');
fprintf('====================================\n');
fprintf('[ OPTIMIZED STORAGE COMPONENTS ]\n');
fprintf('Battery Capacity           : %7.2f kWh\n', E_batt_opt);
fprintf('Fuel Tank Capacity         : %7.2f Liters\n', V_fuel_opt);

fprintf('\n[ FIXED HARDWARE COMPONENTS ]\n');
fprintf('Integrated APU Output      : %7.2f kW\n', data.P_eg_fixed);
fprintf('Traction Motor Output      : %7.2f kW\n', data.P_motor_fixed);

fprintf('\n[ PHYSICAL WEIGHT BREAKDOWN ]\n');
fprintf('Glider/Chassis Mass        : %7.1f kg\n', data.m_glider);
fprintf('Battery Pack Weight        : %7.1f kg\n', mass_batt);
fprintf('Fuel + Tank Weight         : %7.1f kg\n', mass_fuel);
fprintf('Engine + Gen System Weight : %7.1f kg\n', mass_eg);
fprintf('Traction Motor Weight      : %7.1f kg\n', mass_motor);
fprintf('-----------------------------------------------------------------\n');
fprintf('TOTAL CURB WEIGHT          : %7.0f kg\n', opt_mass);

fprintf('\n[ PACKAGING & AREA CONSTRAINTS ]\n');
fprintf('Battery Area               : %7.3f m^2 (Derived via Formula)\n', A_batt_opt);
fprintf('Fuel Tank Area             : %7.3f m^2\n', A_tank_opt);
fprintf('Engine + Generator Area    : %7.3f m^2 (Base: %.3f m^2)\n', A_eg_opt, data.A_eg_base);
fprintf('Traction Motor Area        : %7.3f m^2 (Base: %.3f m^2)\n', A_motor_opt, data.A_motor_base);
fprintf('-----------------------------------------------------------------\n');
fprintf('TOTAL POWERTRAIN AREA      : %7.3f m^2  (Constraint Max: %.2f m^2)\n', opt_Area, data.A_max_limit);

fprintf('\n[ ENERGY CONSUMPTION & RANGE PERFORMANCE ]\n');
fprintf('Calculated EV Consumption  : %7.2f kWh/100km\n', Ec_EV_opt);
fprintf('Pure EV Range (WLTC)       : %7.1f km   (Constraint Min: %.0f km)\n', R_EV_opt, data.min_ev_range);
fprintf('Total Extended Range       : %7.1f km\n', max_range);
fprintf('=================================================================\n\n');

%% =========================================================================
%% SOLVER FUNCTIONS
%% =========================================================================

function [neg_Range, Range] = obj_range(x, data)
    V_fuel = x(1); E_batt = x(2);
    
    m_batt = (data.g_b_slope * E_batt) + data.g_b_fixed;
    m = data.m_glider + (data.g_f*V_fuel) + m_batt + ...
        (data.g_eg*data.P_eg_fixed) + (data.g_m*data.P_motor_fixed);
    
    v_norm = 100;
    P_roll = (m * data.g * data.f * v_norm) / 3600;
    P_req = (P_roll + (data.Cd * data.A_front * v_norm^3) / 76140) / data.eta_t;
    E_c = P_req * 1; 
    
    E_c_EV = E_c + data.P_aux;
    F_c_100km = E_c_EV / (data.rho_Ex * data.eta_eng * data.eta_gen);
    R_fuel = (V_fuel * 100) / F_c_100km;
    R_EV = ((E_batt * data.DoD) / E_c_EV) * 100;
    
    Range = R_EV + R_fuel;
    neg_Range = -Range;
end

function [c, ceq] = dynamic_constraints(x, data)
    V_fuel = x(1); E_batt = x(2);
    
    m_batt = (data.g_b_slope * E_batt) + data.g_b_fixed;
    m = data.m_glider + (data.g_f*V_fuel) + m_batt + ...
        (data.g_eg*data.P_eg_fixed) + (data.g_m*data.P_motor_fixed);
    
    % 1. Footprint Area Calculations
    A_tank  = data.A_tank_base * (V_fuel / data.V_base);
    A_batt  = (E_batt + 1.42) / 23.30; % EMPIRICAL EQUATION APPLIED HERE
    A_eg    = data.A_eg_base   * (data.P_eg_fixed / data.P_eng_base);
    A_motor = data.A_motor_base * (data.P_motor_fixed / data.P_motor_base);
    A_total = A_tank + A_batt + A_eg + A_motor;
    
    % 2. 120 km/h @ 6% Incline Gradeability Performance Check
    v_cruise = 120; incline = 0.06;
    P_cruise = (1/data.eta_t) * ( (m*data.g*data.f*v_cruise)/3600 + (m*data.g*incline*v_cruise)/3600 + (data.Cd*data.A_front*v_cruise^3)/76140 );
    P_eng_req = P_cruise + data.P_aux;
    
    % 3. Dynamic Energy Consumption & EV Range Evaluation
    v_norm = 100;
    P_req_norm = (1/data.eta_t) * ( ((m*data.g*data.f*v_norm)/3600) + ((data.Cd*data.A_front*v_norm^3)/76140) );
    E_c = P_req_norm * 1;
    E_c_EV = E_c + data.P_aux;
    R_EV = ((E_batt * data.DoD) / E_c_EV) * 100;
    
    % Assign inequality constraints: c(x) <= 0
    c(1) = A_total - data.A_max_limit;       % Spatial Footprint Boundary
    c(2) = P_eng_req - data.P_eg_fixed;      % Gradeability Powertrain Sufficiency
    c(3) = data.min_ev_range - R_EV;         % Hard target for pure EV performance
    
    ceq = []; 
end