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?