Finding a suitable algorithm to solve my problem

I want to find an algorithm that can be used to solve the reservoir scheduling problem. My problem has some unbalanced constraints, and it appears to be a non-linear problem. However, I've tried several algorithms and haven't been able to find a suitable one to solve this problem. My code is as follows:

# 构建水库调度的优化函数
using Statistics, Optimization, OptimizationNLopt, Interpolations

# 定义reservoir对象
default_levels = Float32.([22, 24, 26, 28, 30, 32])
default_storages = Float32.([0.61e8, 0.694e8, 0.876e8, 1.133e8, 1.450e8, 1.852e8])
inflow = Float32.([28, 69, 105, 367, 1320, 2440, 2760, 3020, 3140, 2900, 2750, 1870, 1300, 1100, 980, 820])
dt = 1

# 约束条件设置
min_outflow = zeros(Float32, (size(inflow, 1)));
max_outflow = ones(Float32, (size(inflow, 1))) .* 2000.0;

max_levels = 26.0;
min_levels = 22.0;

# 定义水库特征曲线及其插值方法
struct Curve
    levels::Vector{Real}
    storages::Vector{Real}
end

curve = Curve(default_levels, default_storages)

function interp(c::Curve, input_series, series_name)
    interp_linear = Nothing
    if series_name == "levels"
        interp_linear = LinearInterpolation(c.levels, c.storages, extrapolation_bc=Line())
    else
        interp_linear = LinearInterpolation(c.storages, c.levels, extrapolation_bc=Line())
    end
    interp_series = interp_linear.(input_series)
    return interp_series
end

re = interp(curve, [23.0, 24.0], "levels")

function water_balance(inflow, outflow, dt)
    Δstorage = (inflow - outflow) * dt * 3600
    return Δstorage
end

function objective(u, p)
    tmp_storages = water_balance.(inflow, u, dt)
    tmp_levels = interp(curve, tmp_storages, "storages")
    minimum(tmp_levels)
end

function constrain(res, u, p)
    min_storage = minimum(water_balance.(inflow, u, dt))
    max_storage = maximum(water_balance.(inflow, u, dt))
    res .= [interp(curve, max_storage, "storages"), interp(curve, min_storage, "storages")]
end

optprob = OptimizationFunction(objective, cons=constrain)
prob = OptimizationProblem(optprob, max_outflow, [1, 1], lb=min_outflow, ub=max_outflow, lcons=[min_levels, -Inf], ucons=[Inf, max_levels])

using OptimizationNLopt
sol = solve(prob, NLopt.GN_ORIG_DIRECT())  
# ERROR: The algorithm Algorithm does not support constraints. Either remove the `cons` function passed to `OptimizationFunction` or use a different algorithm.

using OptimizationOptimJL
sol = solve(prob, NelderMead()) 
# The algorithm NelderMead{Optim.AffineSimplexer, Optim.AdaptiveParameters} does not support constraints. Either remove the `cons` function passed to `OptimizationFunction` or use a different algorithm.

using OptimizationBBO
sol = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited()) 
# The algorithm BBO_adaptive_de_rand_1_bin_radiuslimited does not support constraints. Either remove the `cons` function passed to `OptimizationFunction` or use a different algorithm.

using OptimizationMOI
optprob = OptimizationFunction(objective, Optimization.AutoForwardDiff(), cons = constrain)
prob = OptimizationProblem(optprob, max_outflow, [1, 1], lb=min_outflow, ub=max_outflow, lcons=[min_levels, -Inf], ucons=[Inf, max_levels])
sol = solve(prob, IPNewton()) 
# ERROR: BoundsError: attempt to access 6-element interpolate((::Vector{Real},), ::Vector{Real}, Gridded(Linear())) with element type Float64 at index [NaN]

SciML / Optimization.jl

Finding a suitable algorithm to solve my problem #610