See PATHSolver.jl/README.md at master · chkwon/PATHSolver.jl · GitHub
But now that I’ve looked at your example a little more deeply, I’m not sure that it is correct, because it doesn’t account for the different i and j at the same x.
How about something like:
function example_one_go(
J::Int,
S::Int,
A::Matrix,
Guess_y::Vector{Float64},
z::Vector{Float64},
C::Matrix,
Guess::Matrix,
)
denom_i, v, last = zeros(S), zeros(S), nothing
function update_cache(i, j, x)
if (i, j, x) != last
denom_i .= A * collect(x[i:J:end])
for s in 1:S
v[s] = A[s, j] / denom_i[s]^2
end
last = (i, j, x)
end
return
end
function f(x_i, x_j, x...)
i, j = round(Int, x_i), round(Int, x_j)
update_cache(i, j, x)
return sum(A[s, j] / denom_i[s] for s in 1:S) / S
end
function ∇f(g, x_i, x_j, x...)
i, j = round(Int, x_i), round(Int, x_j)
update_cache(i, j, x)
fill!(g, 0.0)
for q in 1:J
g[2+J*(q-1)+i] = -sum(A[s, q] * v[s] for s in 1:S) / S
end
return
end
model = Model(PATHSolver.Optimizer)
@operator(model, op_f, 2 + J * J, f, ∇f)
@variable(model, x[i = 1:J, j = 1:J] >= 0, start = Guess[i, j])
@variable(model, y[i = 1:J] >= 0, start = Guess_y[i])
fix(y[1], 1.0; force = true)
@constraints(model, begin
[i in 1:J, j in 1:J], -y[i] * op_f(i, j, x...) + y[j] * C[i, j] / z[j] ⟂ x[i, j]
# 0.0 ⟂ y[1] # Added by default
[i in 2:J], -sum(x[j, i] * C[i, j] for j in 1:J) + z[i] ⟂ y[i]
end)
optimize!(model)
@assert is_solved_and_feasible(model)
X, Y = value.(x), value.(y)
return X, Y, Y ./ z .* C .* X
end