Thanks for the response. Consider the following code:
A = rand(1,10)
B = rand(1,10)
T=length(A)
M=10^6
c_S = 60
c_R = 100
using JuMP
import Juniper
import HiGHS
import MultiObjectiveAlgorithms as MOA
model = JuMP.Model(() -> MOA.Optimizer(HiGHS.Optimizer))
set_attribute(model, MOA.Algorithm(), MOA.Dichotomy())
set_attribute(model, MOA.SolutionLimit(), 5)
@variable(model, S >= 0);
@variable(model, 1 <= R <= 3);
@variable(model, y1[t = 1:T], Bin);
@variable(model, y2[t = 1:T], Bin);
@variable(model, Y[t = 1:T] >= 0);
@variable(model, O[t = 1:T] >= 0);
@variable(model, V[t = 0:T] >= 0);
@variable(model, Q[t = 1:T] >= 0);
@constraint(model, V[0] == 0);
@constraint(model, [t in 1:T], Q[t] == B[t] * R);
@constraint(model, [t in 1:T], Y[t] <= A[t]);
@constraint(model, [t in 1:T], Y[t] <= V[t-1]);
@constraint(model, [t in 1:T], Y[t] >= A[t] - M * (1 - y1[t]));
@constraint(model, [t in 1:T], Y[t] >= V[t-1] - M * y1[t]);
@constraint(model, [t in 1:T], V[t] <= V[t-1] + Q[t] - O[t] - Y[t]);
@constraint(model, [t in 1:T], V[t] <= S - Y[t]);
@constraint(model, [t in 1:T], V[t] >= V[t-1] + Q[t] - O[t] - Y[t] - M * (1 - y2[t]));
@constraint(model, [t in 1:T], V[t] >= S - Y[t] - M * y2[t]);
@objective(model, Max, [-c_S * S - c_R * R, sum(Y[t] for t in 1:T)])
optimize!(model)
solution_summary(model)
println(objective_value(model))
println("S = ", value.(S), "\n","A = ",value.(R))
Now I want to have a solution when I am changing the first objective function to the quadratic term -c_S * S^2 - c_R * R (similar to the approach of @Arpit_Bhatia in NL MO portfolio optimization).