A piece-wise defined objective function can be solved by Gurobi?

How about something like this:

using JuMP
import Ipopt

function MoM(T, R)
    det_R = [r - R[i] for (i, r) in enumerate(R[2:end])]
    det_w = [t - T[j] for (j, t) in enumerate(T[2:end])]
    me_R = (R[end] - R[1]) / T[end]
    A = R[end] * (1 - sum(det_w.^2) / sum(det_w)^2)
    B = sum((det_R - me_R * det_w).^2)
    # Based on the equations from reference, calculate the estimates for alpha 
    # and beta
    beta_est = B / A
    alpha_est = me_R / beta_est
    return alpha_est, beta_est
end

function main()
    T_fc = [0, 185, 348, 515, 658, 830]
    R_fc = [0.1803, 0.1978, 0.2033, 0.2092, 0.2157, 0.2254]
    R0 = R_fc[1]
    Imin, Inom, Imax = 0.2, 0.7, 2
    Lmin, Lnom, Lmax = 0.8305, 2.38, 3.1
    FT = 0.1803 * (1+5.39*0.1)
    Rfc1 = 0.1803 + 0.01
    Rfc2 = 0.1803
    L11, L21 = 2.1, 2.1
    Dindf = [250, 250]
    Ldf = [5, 4]
    C = 20
    α_est, β_est = MoM(T_fc, R_fc)
    α_0, β_0 = α_est / C, β_est * C
    K_dL = 30 * 5.39 # constant to tune load varying casued R aging rate -> KdL
    Dv_dL = K_dL * 5.93e-7 * R0 / (Lmax - Lmin)
    w1 = (FT - Rfc2) / (2FT - Rfc1 - Rfc2)
    w2 = 1 - w1
    a, b = 0.0019727939, 0.0078887
    lb = maximum(Lmin ./ Ldf)
    ub = minimum(Lmax ./ Ldf)
    # Model
    model = Model(Ipopt.Optimizer)
    @variable(model, lb <= x <= ub, start = 0.5)
    @variable(model, lb <= y <= ub, start = 0.5)
    @constraint(model, x + y == 1)
    @NLexpression(
        model,
        f0dL,
        (w1 * (x * Ldf[1] - L11)^2 + w2 * (y * Ldf[1] - L21)^2) * Dv_dL * 3600,
    )
    @NLexpression(
        model,
        α_x,
        ifelse(
            x * Ldf[1] < Lnom,
            a * (x * Ldf[1] - Lnom)^2 + α_0,
            b * (x * Ldf[1] - Lnom)^2 + α_0,
        ),
    )
    @NLexpression(
        model,
        α_y,
        ifelse(
            y * Ldf[1] < Lnom,
            a * (y * Ldf[1] - Lnom)^2 + α_0,
            b * (y * Ldf[1] - Lnom)^2 + α_0,
        ),
    )
    @NLobjective(model, Min, Dindf[1] * (w1 * α_x  + w2 * α_y) + f0dL)
    optimize!(model)
    return value(x*Ldf[1]), value(y*Ldf[1])
end
main()

I think you have a few typos, for example it should probably be β_0 in the ifelse?