Nested optimisation problem in JuMP

No problem.

I’m not sure your code is quite right for the TSP. Here’s how I would write your code. It uses the MTZ formulation of a TSP.

Hopefully it gives you some ideas:

using JuMP
import HiGHS
import MultiObjectiveAlgorithms as MOA
demand_volume = Float64[9,258,33,73,250,8,39,75,80,100]
supply_capacity = [400, 250, 230]
supply_type1_binary = [1, 0, 1]
supply_type2_binary = [0, 1, 0]
number_of_demands = length(demand_volume)
number_of_sources = length(supply_capacity)
distances = [
    0  10 15 20
    10  0 35 25
    15 35  0 30
    20 25 30  0
]
n = size(distances, 1)
model = Model(() -> MOA.Optimizer(HiGHS.Optimizer))
@variables(model, begin
    x[1:number_of_demands, 1:number_of_sources], Bin
    y[1:n, 1:n], Bin
    2 <= u[1:n] <= n
end)
fix(u[1], 1; force = true)
@constraints(model, begin
    [i in 1:number_of_demands], sum(x[i,:]) <= 1
    demand_volume' * x .<= supply_capacity
    [i in 1:n], sum(y[i, :]) == 1
    [i in 1:n], sum(y[:, i]) == 1
    [i in 1:n, j in 1:n; i != 1 && j != 1],
        u[i] - u[j] + 1 <= (n - 1) * (1 - y[i, j])
end)
@expressions(model, begin
    p1_expr, demand_volume' * x * supply_type1_binary
    p2_expr, demand_volume' * x * supply_type2_binary
    p3_expr, sum(distances .* y)
end)
@objective(model, Max, [p1_expr, p2_expr, -p3_expr])
optimize!(model)
solution_summary(model)